STORM 1.0: a simple, ﬂexible, and parsimonious stochastic rainfall generator for simulating climate and climate change

. Assessments of water balance changes, watershed response, and landscape evolution to climate change require representation of spatially and temporally varying rainfall ﬁelds over a drainage basin, as well as the ﬂexibility to simply modify key driving climate variables (evaporative demand, overall wetness, storminess). An empirical–stochastic approach to the problem of rainstorm simulation enables statistical realism and the creation of multiple ensembles that allow for statistical characterization and/or time series of the driving rainfall over a ﬁne grid for any climate scenario. Here, we provide details on the STOchastic Rainfall Model (STORM), which uses this approach to simulate drainage basin rainfall. STORM simulates individual storms based on Monte Carlo selection from probability density functions (PDFs) of storm area, storm duration, storm intensity at the core, and storm center location. The model accounts for seasonality, orography, and the probability of storm intensity for a given storm duration. STORM also generates time series of potential evapotranspiration (PET), which are required for most physically based applications. We explain how the model works and demonstrate its ability to simulate observed historical rainfall characteristics for a small water-shed in southeast Arizona. We explain the data requirements


Introduction
Models of watershed response (rainfall-runoff), water balance, and landscape evolution require characterization of the driving climate, particularly spatially explicit rainfall fields over a time series. The spatial and temporal variability in water delivery to the land surface from the sky and its fate within a drainage basin are a major control on a. the partitioning of water between infiltration and runoff, which affects flood risk and water resources (Beven and Freer, 2001;Beven et al., 1995;Michaelides and Wilson, 2007;Slater et al., 2015;Liang et al., 1994); appropriate spatial representativeness or length of record to well characterize spatial heterogeneity, although radar data can be helpful to improve spatial representation of storms. Second, rain gauge data represent only one realization of multiple potential temporal sequences and spatial patterns of rainfall. Third, global climate models (general circulation models, GCMs) operate at spatial resolutions that are too coarse to represent heterogeneous rainfall fields over small basins or rainfall intermittency (Grotch and MacCracken, 1991;Trenberth et al., 2017). Fourth, while weather generators (regional circulation models, RCMs, or convection permitting models, CPMs) can downscale GCM output for use in dynamic simulation of weather, the model interaction is unidirectional and a regional model is wholly reliant on the boundary conditions provided by the GCM, so it may not well characterize the regional dynamics of climate change (Prein et al., 2015(Prein et al., , 2017Endris et al., 2013;Dunning et al., 2017). Fifth, weather generators that operate at high spatial resolution based on the relevant physics of air mass movement and precipitation formation require detailed information on winds and storm trajectories (Skamarock et al., 2008) that are not available for most basins, and which are challenging to summarize over longer periods of time (e.g., decades). A further consideration is that internationally agreed climate change scenarios themselves (e.g., the Intergovernmental Panel on Climate Change, IPCC) may limit our learning about the regional expression of recent and future climate change because they have already constrained the problem into covarying sets of gridded climate variables from GCM output (or reanalysis data products) based on particular global emissions scenarios. Given major challenges of GCMs to represent rainfall at the basin scale with variance in topography and orography, ensemble output of climate change projections from GCMs (e.g., CMIP5, CMIP6) is unlikely to provide good characterization of the regional or local expression of climate change (Dunning et al., 2017). More importantly, GCM output and reanalysis data products do not allow for the flexibility to assess watershed responses to a wide range of potential regional climate changes that could impact runoff/flood regimes, groundwater recharge, the water balance between plants and the hydrologic cycle, and basin-wide erosion and topographic development. This is especially true in regions where orography and other complicating land surface dynamics affect rainfall fields.
There is an existing class of rainfall space-time rainfall generator models, which are capable of creating spatially explicit rainfall fields for various purposes (Paschalis et al., 2013;Peleg and Morin, 2014;Niemi Tero et al., 2015;Peleg et al., 2017;Benoit et al., 2018). However, each has its peculiarities in operation, data requirements, computational efficiency, programming language, and output resolution, limiting general applications of rainfall simulation to a wider range of modeling applications such as watershed response, ecohydrology, geomorphic landscape evolution, and land-surface responses to a changing climate.
To fill this research gap, we have developed the STOchastic Rainfall Model (STORM), which generates time series of spatially explicit rainfall fields over a gridded domain and spatially uniform time series of evaporative demand. Together, these time series can be used to drive watershed responses within dynamic rainfall-runoff models, water balances within land surface models, or runoff/erosion regimes within landscape evolution models. STORM was introduced elsewhere to explore rainfall patterns and processes within a small drainage basin in the southwest US (Singer and Michaelides, 2017), but here we provide the relevant detail about the model initialization, operation, and evaluation, and describe various improvements that have been made since its initial appearance in the scientific literature. Additionally, we suggest several additional modifications to STORM that may improve its long-term utility, and we outline its potential for broader use in various hydrologic, ecohydrologic, and geomorphic applications. We also provide links to open-source code for STORM in two forms (Matlab and Python), along with sample input data and parameters.

Model initialization and operation
STORM is an empirical-stochastic rainfall generator designed for simple, heuristic simulation of high-resolution drainage basin rainfall under control climate conditions or under different classes of climate change. The term "empirical-stochastic" (sensu, Singer and Dunne, 2004) refers to Monte Carlo selection of several key rainstorm characteristics from distributions that are created from historic datasets. STORM performs this multilayer parameter selection to create multiple sequences of spatially varying rainfall over a drainage basin and over a multidecadal time series. STORM output is particularly useful for simulating rainfall patterns in regions subjected to convective rainfall, where gridded datasets of precipitation do not capture the dynamic behavior (spatial and temporal variability) of rainstorms, and where GCM output is of limited utility in understanding expected changes in rainfall regimes over a basin.
Implementing STORM involves initialization and operation steps (Fig. 1). Initialization includes the creation of distributions of relevant variables that characterize rainstorms and the rainfall they produce over a spatial grid. STORM can be run with one or more hydrologic seasons with different rainfall characteristics. Model inputs include the following (Fig. 2): annual/seasonal precipitation total, P Total (Fig. 2a); storm area, which determines which grid locations in the basin are "hit" by each storm (Fig. 2b); storm center location on a storm grid with defined spacings (Fig. 2c); storm duration, P D (Fig. 2d); rainfall intensity-duration curve number based on a family of curves from which intensity, P I , at the storm center is determined based on the selected value of P D (Fig. 2e); and storm spatial gradient, or the decline in P I with distance from the storm center (Fig. 2f). An additional distri-  Figure 1. Schematic flow diagram illustrating key steps in STORM initialization and operation. Initialization involves gathering data, creating input PDFs, generating input parameters, creating and input/output (I/O) folder structure, and selecting a scenario of climate change, as well as the number of simulations and the length of each simulation. Operation proceeds with the selection of a seasonal or annual precipitation total threshold, followed by generation of rainstorms until the total threshold is satisfied, at which point the season or year changes. Finally, output data are written matrices and files. These steps are outlined in detail within the source code documentation: https://github.com/ blissville71/STORM (last access: 6 September 2018).
bution that is required to create a complete dataset for driving other models is potential evapotranspiration (PET). The PET distribution contains historic daytime and nighttime values for the region of interest organized by month of the year. This collection of values is sampled on a twice-daily basis to characterize average daytime and nighttime evaporative demand that follows the same time signature as the gridded rainfall. PET is assumed to be spatially uniform across the watershed for any 12 h period. Each of these initializing variable distributions can be created from data from a drainage basin, should these data be available (e.g., from one, several, or a network of rain gauges and meteorological stations), or theoretically, based on nearby basin data or generalizing assumptions. Finally, characterizing interarrival times between storms is necessary to enable explicit watershed responses to rainfall inputs (Fig. 3). A distribution of interarrival times is assembled from historic data at all relevant gauging stations. The probability density functions (PDFs) that represent these historical data for each model variable reside external to the STORM code which calls them, and they can thus be easily modified as new information becomes available. For this paper, PDFs were fit manually using Matlab's distribution fitting tool (distfittool), but we recommend that this be automated using a code that optimizes the fit based on maximum likelihood estimators. Intensity-duration curves were generated (annually or for each season) by first computing the maximum P I value for each value of P D (per minute). We then smoothed the resulting points by locally weighted scatterplot smoothing (LOWESS; Cleveland, 1979) using the "smooth" function in Matlab with a span of 0.1 with a second-order polynomial model. Subsequently, we fit a distribution to the LOWESS curve using the curve fitting tool in Matlab. This resulted in the functional form of a double negative exponential curve, consistent with other datasets from dryland regions (Nicholson, 2011). The general equation for the P I −P D curves in Fig 2e is P I = λ×exp(−0.508×P D )+κ ×exp(−0.008×P D )+c. Parameter values for curve 1 are λ = 642.2; κ = 93.1; c = 4.5. Next, we maintained the same functional form of the fitted maxima curve (curve 1) but decreased the magnitude of its coefficients (not the decay parameters) and the intercept, where relevant, by percentiles (90th, 80th, 70th, . . . 10th, 5th) to generate multiple curves that occupy the full phase space of measured P I − P D pairs ( The phase space of P I −P D is not uniform in terms of probability of occurrence of each rainfall event. There is a tendency for more intense storms to occur less frequently than less intense ones. Therefore, each P I − P D curve is assigned a probability of selection that reflects this fact, wherein the most intense curves (1-3) are assigned lower probability of selection (−30 %, −20 %, −10 %, respectively) and the least intense ones (9-11) are assigned higher probability of selection (+30 %, +20 %, +10 %, respectively) with respect to a hypothetical uniform distribution in curve selection probability (Fig. 4a). In other words, rather than assigning a uniform selection probability of 0.091 to all 11 curves, we decreased the probability of selection of the top three and increased the probability of selection of the bottom three curves. This yields the adjusted curve of probabilities labeled as control in Fig. 4a. The empirical curve-fitting method described here is not very elegant or generalizable to different areas, so in future versions of STORM, we will explore the use of cupolas to represent intensity-duration relationships as marginal probabilities (e.g., Vandenberghe et al., 2011).
STORM also incorporates orography, or the tendency for rainfall to vary with elevation, which is common in many drainage basins around the world. We simulate orography by further modifying the probability of curve number selection from the control curve to account for higher intensity rainfall at higher elevations (Fig. 4b). Based on a hypsometric analysis of rainfall as a function of elevation from gauging records, we divided all rainfall grid locations into three orographic groups based on elevation ( For each year, a selection is made from the PDF of P Total (a). Then, during each simulation year, selections are made from the following PDFs on a storm-by-storm basis: storm area (b), storm center location on a Universal Transverse Mercator (UTM) grid, red dots, within a 5 km buffer around the watershed boundary (c), storm P D (d), intensityduration (P I − P D ) curve number (e) based on probabilistic selection favoring less intense storms (see below), and storm intensity gradient with distance from storm center (f).
Singer and Michaelides, 2017). Then we modified P I − P D curve selection probabilities for each storm center accordingly. The lowest orographic group (OG1) has the selection probability of curve 1 (most intense curve) decreased by 50 %, while the mid-elevation group (OG2) has a 25 % reduction of probability for curves 1 and 11, and the highest orographic group (OG3) has a 50 % reduction of probability for selecting curve 11 (least intense curve) (Fig. 4b). This simple procedure of increasing/decreasing the probability of storm intensity at the storm center location appears to capture the general form of orography in the test basin ( Fig. 2d from Singer and Michaelides, 2017). However, a more explicit or theoretical method for characterizing the effects of orographic precipitation could replace the current method in STORM (e.g., including wind speed and direction, as relevant). We imagine that the cupola method mentioned above would also be suitable for characterizing orography by fitting separate cupolas to P I −P D data from gauges within different bands of elevation in a basin determined by hypsometry. STORM is implemented as a function in Matlab with the following syntax:  The input arguments are as follows. MODE refers to whether STORM is running in validation or simulation model. NUMSIMS refers to the number of n-year Monte Carlo simulations to be run. NUMSIMYRS is the number of years in each simulation. SEASONS indicates either one or two seasons with different PDFs to sample from. PTOT_SCENARIO refers to the climate change scenario to be simulated with respect to total annual or seasonal rainfall (wetness) in season 1. STORMINESS_SCENARIO refers to the climate change scenario to be simulated with respect to rainfall intensity (storminess) in season 1. PTOT_SCENARIO2 refers to the climate change scenario to be simulated with respect to total annual or seasonal rainfall (wetness) in season 2. STORMINESS_SCENARIO2 refers to the climate change scenario to be simulated with respect to rainfall intensity (storminess) in season 2. ET_SCENARIO refers to climate change scenario for evapotranspiration (evaporative demand). Each of these climate change scenarios for rainfall can be implemented as either step changes (up or down) or as temporal trends (up or down) playing out over multiple decades, and there is full flexibility to modify the magnitude of these changes for both seasons. Currently, the ET climate change scenario only permits step changes. We describe the climate change scenarios in more detail below. Current climate conditions can also be simulated for any or all of these input parameters.
In Python, STORM is implemented similarly but with some syntax differences. It is likewise implemented as a function, which can be simply imported from the defining script. This function is defined as storm(mode, numsims, numsimyrs, seasons, ptot_scenario, storminess_scenario, ptot_scenario2, storminess_scenario2, ET_scenario, storminess_scaling_factor = 0.05, storm_stepchange = 0.25, storminess_scaling_factor2 = 0.05, storm_stepchange2 = 0.25, ±30 %, ±20 %, and ±10 % probability for the first/last three curves to reflect the fact that larger storms are less probable than smaller ones (a). This initial case is modified as follows to account for basin orography (b). To represent orographic group 1 (OG1 -lowest elevation group, Fig. 2c from Singer and Michaelides, 2017), we decreased the probability of curve 1 by 50 % and applied the difference uniformly to all other curves such that total probability equals 1. To represent orographic group 2 (OG2 -middle elevation group), we decreased the probability of curves 1 and 11 by 25 % and applied the difference uniformly to all other curves. To represent orographic group 3 (OG3 -highest elevation group), we decreased the probability of curve 11 by 50 % and applied the probability difference uniformly to all other curves. Intensity-duration curve numbers correspond to those listed within circles in Fig. 2e. ptot_scaling_factor = 0.05, ptot_scaling_factor2 = 0.05, ET_scaling_factor = 0.25).
The first nine arguments are identical to their Matlab equivalents. The remaining arguments that are supplied with default values permit direct control through the function call of the step and gradual change values specified by the preceding scenarios (see below); in the Matlab version, these are controlled directly from the script defining the function. Documentation is provided as a docstring to the function. The Python version of STORM has a number of very common Python packages as dependencies -os, time, numpy, datetime, matplotlib, scipy, and six. These are typically installed as standard in Python distributions. STORM in Python also requires the pyshp package, which is readily available through package managers including pip and conda. STORM in Python is agnostic between Python 2 and Python 3. In Python, input files (see below) are supplied as comma separated value (csv) files, rather than .mat files, but the contents of these files are identical to those in Matlab. Where distribution objects are provided in Matlab, these are provided to the Python version as simple text files specifying the key parameters of those distributions. More details of both of these input file formats are provided in the docstrings.
To operate STORM, we simulate rainfall on a storm-bystorm basis with a temporal resolution of 1 min at each rainfall output grid location with a resolution of 1 km. The high spatial and temporal resolution enables rich information to be generated from STORM, allowing for detailed assessment of spatial heterogeneity and temporal variance in rainfall fields. We first select a threshold value of P Total for a simulation year from its distribution and then generate rainstorms until the median running total of P Total across all rainfall output grid locations within the basin equals or exceeds the selected threshold value of P Total . Then a new simulation year begins and STORM proceeds until the length of simulation (e.g., several decades) is complete. P Total can be designated as either a seasonal or annual total, depending on whether the watershed of interest is characterized by strong seasonality in rainfall. If STORM is implemented with two seasons, separate PDFs of rainstorm characteristics (Fig. 2a, b, d, and e) should be prepared for each season. PDFs for this paper were fit to historical data from the Walnut Gulch Experimental Watershed (WGEW) using the distribution fitting tool within Matlab v2017b. Sample distributions are provided at https://github.com/blissville71/STORM (last access: 6 September 2018).
To model climate change as a step change in wetness (shift in total annual precipitation), we shifted the P Total distribution up P Total (+) and down P Total (−) by 1 standard deviation without changing its shape (Fig 2a). To model climate change as a step change in storminess (shift in intensity for a particular rainstorm duration), we modified the selected intensity for all storms as P I ± × P I , where is a fractional change in storm intensity (Fig 2e). Thus, the selected intensity for all rainstorms is raised for P I (+) and lowered for P I (−). We have also built into STORM the capability to assess trends in both P Total and P I by multiplying the selected value of either rainfall variable for each year of simulation by an annual change scalar (see below for details).
The results generated at each rainfall output grid location over all simulations can be statistically analyzed. STORM generates a rainstorm matrix that includes storm no., storm area (km 2 ), storm duration (min), intensity-duration curve no., storm intensity at storm center (mm h −1 ), no. of gauging (or rainfall grid) locations hit, intensity recession value (mm h −1 km −1 ), storm total (mm), longitude (m), latitude (m), year (yr), and cumulative simulation time (h). STORM also generates a separate matrix of output at each rainfall grid location that includes year, storm no., local storm intensity (mm h −1 ), storm duration (min), local storm total (mm), annual local cumulative precipitation total (mm), interarrival time between storms (h), and cumulative simulation time (h).
All simulations of STORM for this paper were done using Matlab v.2017b on an Intel Dual Core i7-4712HQ CPU at 2.30 GHz running on a Dell Inspiron laptop with 16 GB of RAM. In simulation mode, each simulation of 30 years on a grid of 128 rainfall locations within a 149 km 2 basin took ∼ 12.5 min for the control run, or ∼ 25 s per simulated year. The computational time of each ensemble varies depending on the scenario. For example, the simulation of P Total (+)P I (−) took ∼ 20 min per 30-year simulation (or ∼ 40 s per simulated year).

STORM application to Walnut Gulch
Here, we provide a few more details about how we assembled the relevant PDFs for WGEW that are shown in Fig. 2. PDFs were generated within Matlab using the distribution fitting tool, though the methods are straightforward enough to be easily applied within many data analysis software packages. To create the PDF of storm areas, we truncated the distribution of Syed et al. (2003) to exclude all rainstorms in their largest spatial bin (which comprised approximately one-half of all their rainstorms) and then fitted an extreme value distribution to the remaining storms (Fig. 2b). Syed et al. (2003) indicated that this bin of high frequency for large spatial ar- eas is probably an artifact of sampling such that they are overestimating storm areas due to repeat sampling of moving storms and due to the occurrence multiple simultaneous storm cells. Thus, our exclusion of this largest area bin is justified to achieve our goal of modeling individual (discrete) rainstorms and the resulting distribution has a mean storm area of just under 90 km 2 . In our initial tests of STORM, we found that the results were very sensitive to the β parameter (Eq. 1; Fig. 2f). Thus, to enable a variety of plausible storm gradients, we allowed β to vary according to a normal distribution (µ = 0.25, σ = 0.08, truncated to the interval [0.15, 0.67 km −1 ]), based on the range of values reported elsewhere for WGEW (Eagleson et al., 1987;Morin et al., 2005).
We analyzed all the rainfall data from WGEW to assess the frequency of different storm intensities and found that, as expected, there is a high frequency of rainstorms with low intensity, and vice versa (Fig. 5). This validates our development of the control set of curve selection probabilities (Fig. 4a).
In each multiyear simulation, storms of various areas cover different parts of STORM's model domain. If we aggregate all the simulated rainfall data from each of the rainfall output grid locations, we see increasing coverage of the P I − P D phase space with increasing length of simulation in years (Fig. 6). Simulations of 25 years appear to fill in much of the P I − P D phase space for control (no-climate-change) conditions (Fig. 6f). This indicates that STORM is faithfully representing the input data on intensity-duration of rainfall in a stochastic treatment over multiple decades.
Finally, a plausible characterization of climate change is necessary to gain insight into the potential impacts to watershed response. In terms of precipitation, climate change can manifest in wetter or drier conditions (e.g., over a season or a year) and/or in different storm characteristics (e.g., relationship between P I and P D ). We therefore suggest there are two classes of climate change that affect convective precipitation and watershed response. In STORM, we can simulate different classes of climate change based on very simple rules. To simulate step changes in basin wetness (annual/seasonal precipitation totals), we shift the PDF of P Total up for P Total (+) or down for P Total (−) by 1 standard deviation, while retaining the same shape of the distribution (Fig. 2a). Note: this truncates the left tail of original P Total distribution for the scenario of P Total (−). To simulate step changes in storminess, we multiply the selected value of P I at the storm center by a scalar fraction and add (subtract) the product to the selected P I value to reflect increased (decreased) storminess (Fig. 2e). Specifically, we modified the selected intensity for all storms as P I ± ×P I , where is a fractional step change in storm intensity. Thus, the selected intensity for all rainstorms is raised for P I (+) and lowered for P I (−). In this paper, we used a step change of = 0.25.
We also enable the simulation of temporal trends in rainfall. We include scalar multipliers that can be set within the STORM code, and which are used to modify the relevant value of rainfall each year. For example, to characterize a trend in wetness (total precipitation per season or year), we annually update the P Total PDF by generating a new PDF for each year after progressively modifying the mean value of the PDF (Fig. 2a) as µ = µ ± µ × , where is a fractional scaling factor that increases the PDF of P Total each year of simulation. Within STORM, a similar procedure is used for trends in storminess. In this case, to characterize a positive trend in storminess (magnitude of storm intensity), we annually update selected P I values by multiplying by a scalar as P I ± × P I , where is a fractional scaling trend that increases the selected value of storm intensity by an accumulating trend each year of simulation (i.e., the storminess trend in any year of simulation is computed as Y = ×Y , where Y is the simulation year). In this paper, we used initial values for both and of 0.05. We note that it is possible to simulate trends in storminess separate from trends in wetness, or in combination.

Model output and evaluation
It is possible to extract from STORM detailed output of rainfall characteristics for discrete rainstorms at each rainfall grid location, including storm-by-storm intensity, duration, and storm totals (Fig. 7). Thus, one can develop a localized time series of these rainstorm characteristics for different locations or sub-basins within a drainage basin. Once the simulated interarrival times between each pair of rainstorms (Fig. 3) are added back to the time series and used alongside the simulated time series of PET (Fig. 8), we obtain a temporally explicit climate driver that can be used to drive models of watershed response, land surface models, or to make localized water balance calculations. Figure 9 illustrates output at a single gauging location of explicit time series of storm rainfall for different climatic scenarios (see below). Note that the time series shown in Figs. 7-9 only illustrate output for a single multidecadal simulation. The stochastic capability of STORM allows for generation of multiple n-year simulations. Thus, the STORM output for all such simulations should ideally be analyzed statistically at each representative location to produce representative time series of rainfall variables. Otherwise, they can be used as ensemble inputs to other model frameworks. We evaluated the model's skill at simulating observed rainfall characteristics at 85 gauges in WGEW. We selected three representative variables: number of storms per year, average storm total, and total annual precipitation. We compared observed versus simulated values of these variables at each gauging location in WGEW. Figure 10 shows these as simulated versus observed on 1 : 1 plots, and relevant statistics on model skill are provided. Generally, STORM demonstrates a high level of skill at simulating rainfall characteristics across the domain, without any model tuning. The spread of points around the 1 : 1 line arises from the ensemble of multidecadal simulations, and it will therefore vary from ensemble to ensemble. Generally, it is desirable to have a range of values that span the 1 : 1 line (above and below it). This ensures the model is generating a field of rainstorm characteristics beyond those which exist in the historical record, a notable strength of the Monte Carlo procedure implemented in STORM.
Note: the random selection of model parameters (Figs. 2 and 3) was done using the random function in Matlab (v.2016b), when sampling from a PDF (e.g., P Total ), and using the datasample function when sampling from a vector of data points (storm center location). The random number generator seed was shuffled prior at the beginning of each Matlab session.  Figure 9. Illustration of simulated rainfall at one gauging location for a 25-year simulation with two seasons. The panels are (a) control climate conditions (CC); (b) a step-change increase in P Total for season 1 and a step-change decrease in P Total for season 2 (S 1 (+) / S 2 (−)); (c) a step-change decrease in P Total for season 1 and a step-change increase in P Total for season 2 (S 1 (−) / S 2 (+); (d) positive trends in P Total for both seasons (S 1 (T+) / S 2 (T+); (e) a step-change increase in P I for season 1 and a step-change decrease in P I for season 2; (f) a step-change decrease in P I for season 1 and a step-change increase in P I for season 2; (g) positive trends in P I for both seasons.

STORM in the context of climate change
Regional gridded datasets provide a picture of trends in certain climate variables that are relevant to rainstorms. Figure 11 shows monthly output at WGEW for the 0.5 • CRU TS3.24.01 dataset for temperature anomalies (Fig. 11a) and precipitation anomalies (Fig. 11b). These datasets show a significant increase in temperature of ∼ 2 • C over recent decades (which constitutes the period of the WGEW rainfall record), but there is no clear trend in the monthly precipitation from this dataset. The lack of a trend in monthly precipitation data contrasts with the high-resolution data from WGEW, which showed a long-term increase in P Total but with declining P I (Singer and Michaelides, 2017). Such changes to the hydrology of a basin could have major implications for its runoff regime, water balance, and landscape evolution. This suggests that gridded global datasets are not adequate for investigations of convective rainfall in dryland (and potentially many other) basins. This observation lends support for both detailed data analysis from dryland datasets such as WGEW but also for modeling approaches such as STORM for exploring the potential impacts of climate change on rainfall and watershed responses. Values of RMSE and PBIAS are also shown. No model tuning was performed to achieve these results. We note that observed values of storm total and P Total are slightly overpredicted by our simulations (over these 30 test ensembles, each of 30 years). Nevertheless, these results demonstrate a high level of model skill in reproducing a range of rainfall characteristics that span observed values without any model tuning. Figure 8 shows STORM output of ET showing clear daynight variations, as well as seasonality (monthly day-night ET were sampled from historical data compiled by month). Figure 8a shows the control climate scenario, Fig. 8b shows a step-change increase in ET, and Fig. 8c shows a step-change decrease. Figure 9 illustrates STORM rainfall output at a single gauging point in the basin as complete time series that include interarrival times, so the data can be used to drive  Figure 11. Monthly mean temperature (a) and monthly mean precipitation anomalies (b) for the closest grid location to Walnut Gulch from the 0.5 • Climate Research Unit (CRU) TS3.24.01 dataset (Harris and Jones, 2017). These plots show a recent increase in mean temperature that could influence rainfall P I − P D relationships, irrespective of the lack of long-term trend in precipitation totals from the gridded dataset. Curves are based on LOWESS smoothing (Cleveland, 1979) via the smooth function in Matlab (v.2016b) using a 30 % data span. other models. The figure also shows these time series of rainstorm totals for various climate change scenarios including control climate (Fig. 9a), a step-change increase in total season 1 rainfall with a step-change decrease in total season 2 rainfall (Fig. 9b), a step-change decrease in total season 1 rainfall with a step-change increase in total season 2 rainfall (Fig. 9c), positive trends in total rainfall for both seasons (Fig. 9d), a step-change increase in season 1 rainfall intensity with a step-change decrease in season 2 rainfall intensity (Fig. 9e), a step-change decrease in season 1 rainfall intensity with a step-change increase in season 2 rainfall intensity (Fig. 9f), and positive trends in rainfall intensity for both seasons (Fig. 9g). These outputs illustrate the range of capability included in STORM for simulating an array of different regional expressions of climate change that could have important implications for watershed response. In general, increases in P Total tend to densify the number of rainstorms in the time series (and vice versa), while increases in storminess tend to create more peaked rainstorm totals (cf. Fig. 9d, and g).

STORM data requirements
An important consideration of any model is the data requirements. Obviously, the lower amount of data required enables more widespread model use to tackle environmental problems. However, on the other end of this spectrum, insufficient data can lead to poor model skill. Thus, we aimed to strike a balance that would enable widespread use of STORM with limited data. In terms of storm event rainfall data, we inves- Figure 12. Illustration of the limited sensitivity of STORM to the quantity of input data. Plots of P I versus P D for various percentages of the complete dataset (indicated at the upper right of each subplot) show that even when only 5 % of the original gauge data are included, the P I − P D phase space for the complete dataset is still broadly defined. The number of data points declines across these subpanels from 185 109 (100 %) to 1851 (1 %). These factors indicate that STORM could be reasonably applied in a basin with dramatically less available gauging data, if other storm characteristics can be constrained by other means. tigated how much rain gauge data are required by STORM to well characterize the historical phase space of rainfall characteristics. For WGEW, we plotted event rainfall intensity versus duration for 100 % of the available data, 50 %, 25 %, 10 %, 5 %, and 1 % (Fig. 12). This analysis shows that even at 1 % of the available data (n = 1851 rainfall events), the P I − P D phase space is still clearly delineated, allowing for development of the necessary P I − P D curves (Fig. 2e). Therefore, one could use data from a basin with only one or several rain gauges that have collected event rainfall for a few decades.
Another important data requirement in STORM is a PDF of storm areas (Fig. 2b). This may be a more challenging PDF to develop due to limited data. Options here include analyzing spatial statistics on discrete storms from a network of rain gauges (e.g., Syed et al., 2003), analyzing storm characteristics from rainfall radar images (e.g., Peleg and Morin, 2012), or developing a hypothetical distribution based on regional understanding of mapped rainstorms (e.g., from hyetographs). Finally, the density of the gauging network could have important influence on the storm intensity gradient with distance from the storm center, so the parameters of this relationship may be less certain for less dense gauge spacings.

Extension/modification of STORM for additional applications
For many applications, a complete account of the water balance between precipitation, evapotranspiration, and infiltration is required. Specifically, closing the local water balance requires quantifying the evaporative demand and the length of inter-storm periods, which enable drainage and drying of soil layers. These factors affect the watershed response to subsequent rainfall events. In order to characterize these aspects, we have modified STORM from its original capability (Singer and Michaelides, 2017) to include several new features. First, we have added a PDF of inter-storm periods that is sampled randomly after each storm event. The addition of these inter-storm periods changes STORM output into time series that reflect real time at the Earth's surface (e.g., Fig. 9). Second, we have assembled a PDF of potential evapotranspiration based on measurements of temperature and relative humidity, metrics which are readily available over multiple scales (Fig. 8). Third, we implemented seasonality in rainfall to enable simulations over a single season or year, or over two seasons with distinct differences in precipitation characteristics (distinct PDFs of rainfall in summer compared to winter). The WGEW example implemented for illustration of the model here is one such basin with a strong monsoon season that produces a high percentage of the annual rain and most of the runoff, compared with the winter season dominated by weak frontal storms. These improvements to STORM now make it suitable as a climate driver of other watershed response models that simulate hydrology between slopes and channels (surface runoff, infiltration, streamflow) Wainwright, 2002, 2008;Michaelides and Wilson, 2007), groundwater recharge during and after rainfall events (Beven and Freer, 2001), and interactions between streamflow and alluvial aquifers (Evans et al., 2018). It also enables STORM to be useful in water balance models (e.g., land surface models) to assess water availability to plants through dynamic ecohydrological simulation of plant-climate interactions and water utilization (D'Odorico et al., 2007;Caylor et al., 2006;Laio et al., 2006), as well as energy/carbon fluxes between the land surface and the atmosphere (Best et al., 2011;Bonan, 1996). Finally, STORM can also be used to drive geomorphic models that characterize erosion and deposition processes within drainage basins in response to sequences of rainfall and runoff (Michaelides et al., 2009Michaelides and Martin, 2012;Michaelides and Singer, 2014), and even landscape evolution models that simulate landform development over longer timescales (Tucker and Hancock, 2010;Hobley et al., 2017). Coupling STORM to such models would enable a wide range of interdisciplinary scientists to investigate key problems in the environment that have their origin in the climate system. These range from which water sources are used by plants (Sargeant and Singer, 2016;Evaristo et al., 2015;Evaristo and McDonnell, 2017;Dawson and Ehleringer, 1991) to what is the dominant source and timing of groundwater recharge (Cuthbert et al., 2016;Wheater et al., 2010;Scanlon et al., 2006) to the role of climate in shaping landscape morphology Tucker and Bras, 2000;Tucker and Slingerland, 1997;Michaelides et al., 2018). A version of STORM is under active development in the modular opensource surface process modeling framework Landlab (Hobley et al., 2017), in part to facilitate such future work. Another key area of future work would be to investigate how temporal resolution of rainfall data affects the signal of observed trends in rainfall (e.g., Barbero et al., 2017) and how these might yield different watershed responses.