A Nonstationary Stochastic Rainfall Generator Conditioned on Global Climate Models for Design Flood Analyses in the Mississippi and Other Large River Basins

Existing stochastic rainfall generators (SRGs) are typically limited to relatively small domains due to spatial stationarity assumptions, hindering their usefulness for flood studies in large basins. This study proposes StormLab, an SRG that simulates precipitation events at 6‐hr and 0.03° resolution in the Mississippi River Basin (MRB). The model focuses on winter and spring storms caused by water vapor transport from the Gulf of Mexico—the key flood‐generating storm type in the basin. The model generates anisotropic spatiotemporal noise fields that replicate local precipitation structures from observed data. The noise is transformed into precipitation through parametric distributions conditioned on large‐scale atmospheric fields from a climate model, reflecting spatial and temporal nonstationarity. StormLab can produce multiple realizations that reflect the uncertainty in fine‐scale precipitation arising from a specific large‐scale atmospheric environment. Model parameters were fitted monthly from December–May, based on storms identified from 1979 to 2021 ERA5 reanalysis data and Analysis of Record for Calibration (AORC) precipitation. StormLab then generated 1,000 synthetic years of precipitation events based on 10 CESM2 ensemble simulations. Empirical return levels of simulated annual maxima agree well with AORC data and show an overall increase in 1‐ to 500‐year events in the future period (2022–2050). To our knowledge, this is the first SRG simulating nonstationary, anisotropic high‐resolution precipitation over continental‐scale river basins, demonstrating the value of conditioning such stochastic models on large‐scale atmospheric variables. StormLab provides a wide range of extreme precipitation scenarios for design floods in the MRB and can be further extended to other large river basins.


Introduction
The Mississippi River Basin (MRB) has experienced multiple extreme floods throughout history, including the Spring 1927 Flood (J. A. Smith & Baeck, 2015), the Great Flood of 1993 (Dirmeyer & Brubaker, 1999;Dirmeyer & Kinter, 2009), and the June 2008 Midwest Flood (Budikova et al., 2010;Dirmeyer & Kinter, 2009).These floods were each caused by sequences of severe large-scale rainstorms, all associated with strong water vapor transport from the Gulf of Mexico (Lavers & Villarini, 2013;Moore et al., 2012;Su et al., 2023).The project design flood was developed to estimate the maximum possible peak discharge and water levels in the lower MRB, providing guidelines for flood protection infrastructure and management (Gaines et al., 2019).Historical records of extreme rainstorms were placed in sequences, with modest perturbations to their timing, location, and severity to create hypothetical worst-case rainstorm scenarios (V. A. Myers, 1959).Of these, Hypo-Flood 58A, consisting of the rainstorm in January 1937 followed shortly by rainstorms in January 1950 and then February 1938, produced the largest discharge on the lower Mississippi River and has been used as the design flood since 1955 (Lewis et al., 2019;McWilliams & Hayes, 2017).However, this approach is not without limitations, including using a limited rainstorm record (i.e., pre-1950s) and providing only one worst-case scenario, rather than a range of possible outcomes with occurrence probabilities.With growing evidence that rainstorms are intensifying under future climate change (e.g., Ban et al., 2015;Fischer & Knutti, 2016;Fowler et al., 2021;Matte et al., 2021), historical design floods may fail to reflect the range of possible scenarios and may underestimate future risks (Wright et al., 2014;Yu et al., 2020).
In principle, stochastic rainfall generators (SRGs) can generate synthetic precipitation data that extend short observation records and provide a wide range of extreme scenarios to support flood risk management (Ailliot et al., 2015;Benoit et al., 2018;Wilks & Wilby, 1999).SRGs can be categorized into three main types based on their simulation dimensions: point or area-averaged models, multi-site models, and space-time models.We will focus on the last type, which generates spatiotemporal precipitation data over regular grids.Reviews of other model types can be found in, for example, Waymire and Gupta (1981), Mehrotra et al. (2006), andQin (2011).
Existing space-time SRGs can be further classified into three major categories: (a) cluster-based models, (b) multifractal models, and (c) meta-Gaussian models.These categories are briefly reviewed here.Cluster-based models conceptualize the precipitation system as a hierarchical organization in which groups of small rain "cells" (∼10-50 km 2 ) are embedded within large rainband areas (∼10 3 -10 4 km 2 , e.g., Waymire et al., 1984).These models simulate rainband centers using Poisson processes, with elliptical rain cells randomly generated around each rainband following a clustered point process (e.g., Chen et al., 2021).Properties of these rain cells, including intensity, lifetime, area, and orientation are modeled via distributions fitted to observed data (Northrop, 1998).Total precipitation intensity at any location is calculated by summing the contributions from all existing (and possibly overlapping) rain cells.Cluster-based models were originally developed to simulate precipitation time series at single points (i.e., Neyman-Scott and Bartlett-Lewis models, Rodriguez-Iturbe et al., 1997) and later expanded for space-time simulations through work by Le Cam (1961), Waymire et al. (1984), Cowpertwait (1995), and Northrop (1998), and others.Recent improvements include randomizing model parameters (Kaczmarska et al., 2014;Onof & Wang, 2020) and applying non-homogenous Poisson processes for rain cell arrivals (Burton et al., 2010).This modeling approach attempts to mimic the physical organization of precipitation systems and enables analytical derivation for key precipitation properties like moments and spatial correlation (Cowpertwait, 1995;Leonard et al., 2008).However, cluster-based models struggle to fully describe the precipitation statistical structure across scales (Foufoula-Georgiou & Krajewski, 1995) and typically lack dependence modeling between precipitation properties, for example, rain cell duration and intensity, with implications for simulating extremes (Chen et al., 2021;Paschalis et al., 2013).The multi-level model structure also introduces many parameters, increasing model calibration difficulty (Jothityangkoon et al., 2000).
Multifractal models were developed based on empirical evidence that precipitation systems exhibit similar statistical properties across different temporal and spatial scales (i.e., scale invariance, see V. K. Gupta & Waymire, 1993;Menabde et al., 1997;Schertzer & Lovejoy, 1987).A common modeling approach, known as the discrete multiplicative random cascade (MRC), involves recursively subdividing the simulation domain into small squares and multiplying the precipitation intensity in each square by a random weight drawn from a fitted distribution (Rupp et al., 2012;Serinaldi, 2010;Sharma et al., 2007).A similar process is applied along the time dimension using a different dividing factor to account for anisotropy between space and time (Marsan et al., 1996;Veneziano et al., 2006).Other model variations include continuous MRC (Marsan et al., 1996), temporal Markov chain coupled with spatial MRC (Jothityangkoon et al., 2000;Over & Gupta, 1996), and rainfall cascade decomposition (Raut et al., 2019;Seed et al., 2013).Multifractal models aim to reflect precipitation physical structures and require fewer parameters than cluster-based models (Raut et al., 2018).Nevertheless, these models assume that precipitation fields are spatially isotropic and only represent limited types of space-time scaling properties.These assumptions can introduce bias when precipitation does not exhibit perfect invariance over different temporal and spatial scales (Serinaldi, 2010).
Meta-Gaussian models represent precipitation evolution by simulating an underlying 2D Gaussian noise field that attempts to mimic the spatiotemporal correlation structure of actual precipitation (Baxevani & Lennartsson, 2015).This noise field is then transformed into precipitation amounts using a non-linear function, often a parametric precipitation distribution, generating simulations that are statistically consistent with observed precipitation (Mascaro et al., 2023).Meta-Gaussian models have become increasingly popular in recent years (Papalexiou et al., 2021;Sigrist et al., 2012;Vischel et al., 2011), with major differences in how they represent precipitation spatial structures.Approaches include parametric covariance functions (Papalexiou, 2018;Paschalis et al., 2013), turning band methods (Leblois & Creutin, 2013), fast Fourier transforms (FFTs, Nerini et al., 2017), stochastic partial differential equations (Fuglstad et al., 2015), and kernel convolutions (Fouedjio et al., 2016).A key advantage of meta-Gaussian models is their flexibility in modeling precipitation distributions and spatiotemporal dependence structures, which can yield improved extreme value simulations compared to other SRGs.However, they often involve many parameters and high computational demand for large-scale simulations (Papalexiou & Serinaldi, 2020;Paschalis et al., 2013).Modeling temporal autocorrelation also remains challenging due to its dependence on storm movement and on spatial correlation (Storvik et al., 2002).Beyond the three major categories described, there also exist some other space-time rainfall generators such as variogrambased simulations (Schleiss et al., 2014) and deep learning models (e.g., Leinonen et al., 2021;Rampal et al., 2022;C. Zhang et al., 2020).Another relevant branch of models are stochastic weather generators, which simulate time series of weather regimes and atmospheric variables like temperature and precipitation across multiple sites via stochastic methods such as hidden Markov models or k-nearest neighbor sampling (see Ailliot et al., 2015;R. S. Gupta et al., 2023;Najibi et al., 2021;Sparks et al., 2018).
Most existing SRGs assume the precipitation field is second-order stationary, that is, the mean of precipitation is constant across the region, with covariance between two grid cells depending only on their separation distance rather than their absolute locations (Benoit et al., 2018;Fuglstad et al., 2015;Hristopulos, 2020).This limits simulation domains to relatively small areas where precipitation can be considered homogeneous (Papalexiou et al., 2021;Schleiss et al., 2014;Vischel et al., 2011).Such an assumption precludes the usage of such SRGs for simulating extreme precipitation over continental-scale basins like the MRB.Extreme rainstorms in large basins are likely to arise from different types of weather systems (e.g., tropical cyclones, extratropical cyclones, mesoscale convective systems) exhibiting distinct precipitation distributions and correlation structures (Barlow et al., 2019;B. Liu et al., 2020;Schumacher & Rasmussen, 2020).Precipitation distributions also vary spatially, with regions near moisture sources or subject to orographic lifting showing higher means and variances (Hitchens et al., 2013;Marra et al., 2022).Nonstationarity in spatial correlation structures is clearly evident in highresolution precipitation observations, with patterns varying significantly within different parts of large weather systems like extratropical cyclones (Nerini et al., 2017;Niemi et al., 2014).Moreover, observed increases in extreme precipitation due to anthropogenic climate change suggest temporal nonstationarity in precipitation properties, challenging the usage of temporally stationary distributions of precipitation amount (Gori et al., 2022;Pan et al., 2016).Consequently, stationary SRGs will struggle in several key respects to represent the strong spatiotemporal nonstationarity of extreme precipitation within a large basin.Several efforts have been made to model nonstationary precipitation fields (e.g., Fouedjio, 2017;Fouedjio et al., 2016;Fuglstad et al., 2015;Nerini et al., 2017), but to our knowledge, no SRG has been developed to simulate nonstationary precipitation fields for a basin as large as the MRB.
In this study, we propose the Space-Time nOnstationary Rainfall Model for Large Area Basins (StormLab), an SRG that generates 6-hr, 0.03°resolution rainstorms over the MRB.Following the meta-Gaussian framework, the model generates spatiotemporally correlated noise fields and converts them to actual precipitation patterns via precipitation distribution functions.Unlike some other SRGs that simulate precipitation without conditioning on additional atmospheric variables, StormLab performs conditional simulations where precipitation distributions vary across space and time, dependent on large-scale atmospheric variables from global climate models (GCMs).Such conditioning incorporates spatial and temporal nonstationarity of precipitation and realistically represents how local-scale precipitation is embedded within and driven by the larger-scale atmospheric environment.Nonstationarity in the spatial correlation of precipitation is further addressed by replicating local correlation structures from precipitation observations using a "local window" FFT approach.StormLab was used to generate 1,000 synthetic years of winter and spring rainstorms conditioned on GCM outputs.StormLab can alternatively be understood as a stochastic downscaling approach, that is, establishing statistical relationships between coarse GCM data and local precipitation observations.In this study, "rainstorm" is used to denote a heavy precipitation event, though it may contain snow or mixed-phase precipitation.All input data and stochastic simulations consider liquid water equivalent, without distinction of precipitation phase.The remainder of the paper is organized as follows: Section 2 describes the study basin and data sets.Section 3 details the proposed methods.Results are shown in Section 4, followed by discussion in Section 5. A summary and conclusions are provided in Section 6.

Study Site and Data Processing
The MRB, the largest watershed in North America, has a drainage area of 3.2 million km 2 and consists of five major subbasins: the Arkansas-Red, Missouri, Upper Mississippi, Ohio-Tennessee, and Lower Mississippi (Figure 1a).Major flood events occur during winter and spring (December-May), resulting in significant socioeconomic impacts across the basin (Mutel, 2010;M. F. Myers & White, 1993;A. B. Smith & Katz, 2013).Due to the basin's large and complex drainage network, the locations and spatiotemporal patterns of heavy precipitation greatly influence flood responses and damages throughout the basin (Su et al., 2023).
The Analysis of Record for Calibration (AORC) data set provides reference precipitation data at 6-hr, 0.03°r esolution over the MRB from 1979 to 2021.The data set was created by combining hourly NLDAS-2 and Stage IV precipitation data (Fall et al., 2023), which we resampled to 6-hr intervals in this study.ERA5 reanalysis, produced by the European Center for Medium-Range Weather Forecast, provides hourly, 0.25°fields of largescale atmospheric variable fields (Hersbach et al., 2023); we used precipitation, precipitable water, integrated water transport (IVT), and 850 hPa winds over the study region from 1979 to 2021.ERA5 fields were resampled to 6-hr intervals and spatially upscaled to match the resolution of the CESM2 data set (described next), using the Python package "xESMF" (Zhuang et al., 2020).The CESM2 Large-Ensemble Project provides 100 ensemble members of historical (1850-2014) and future projection (2015-2100; SSP3-7.0 scenario) simulations using the CESM2 model (Danabasoglu et al., 2020;Rodgers et al., 2021).We used 10 ensemble members from 1951 to 2050, equivalent to 1,000 years of atmospheric simulations, from which precipitation, precipitable water, IVT, and 850 hPa winds were retrieved at approximately 1°resolution and 6-hr intervals.Both AORC precipitation and ERA5/CESM2 precipitation represent the total accumulated amount over 6 hr.For example, the value at 06:00 UTC represents the total precipitation from 00:00 to 06:00 UTC.In contrast, precipitable water, IVT, and 850 hPa winds from ERA5/CESM2 represent the average amount over 6 hr.For instance, the precipitable water at 06:00 UTC is the average precipitable water from 00:00 to 06:00 UTC.The CESM2 large-scale atmospheric variables were bias-corrected against ERA5 data for each season using the CDF-t method (M.Vrac et al., 2012).The CDF-t method finds a transformation that maps the cumulative distribution function (CDF) of the CESM2 variable in the historical period  to the corresponding ERA5 CDF (Pierce et al., 2015).The same mapping is then applied to the CDF of CESM2 variable in the early  and future (2022-2050) periods.This preserves trends and changes in the CESM2 data while removing bias (see Supporting Information S1).

Rainstorm Identification
Rainstorms associated with strong water vapor transport were identified over the MRB in ERA5 and CESM2.First, we applied a storm tracking method STARCH (Y.Liu & Wright, 2022b) to identify periods of intense water vapor transport over the tracking domain (29-50°N, 79-113°W, Figure 1a).STARCH applies a high threshold (500 kg m 1 s 1 ) to the IVT fields to identify regions of intense water vapor transport at each time step.These regions were then expanded to a lower threshold (250 kg m 1 s 1 ) to capture more IVT areas.The identified IVT objects were tracked across time by computing the overlap ratio between objects at consecutive time steps.If the overlap ratio exceeds 0.2, the objects were considered the same IVT event and tracked through time (Figure 1a, see Text S2 in Supporting Information S1 for more details).These IVT events were considered as indicators for individual rainstorms and associated with concurrent precipitation, precipitable water, and wind speed fields from ERA5/CESM2 data (e.g., Figure 1b).To focus on extreme events that potentially lead to large-scale flooding, we only selected rainstorms with precipitation covering >10% of the basin for more than 24 hr.The starting and ending periods of a rainstorm were also cropped if the precipitation coverage was <10% of the basin.For ERA5 data, we identified 1,088 rainstorms from December-May during 1979-2021 and associated AORC precipitation (Figure 1c).For CESM2, we identified 26,166 rainstorms over the 1,000 years.

Stochastic Rainfall Generator
The Space-Time nOnstationary Rainfall Model for Large Area Basins (StormLab) was developed based on ideas from Notaro et al. (2014) under a meta-Gaussian framework.For a given ERA5/CESM2 rainstorm, 2D Gaussian noise was generated to represent precipitation spatiotemporal structures and converted to precipitation amounts through parametric distributions conditioned on large-scale atmospheric variables.Multiple noise realizations were simulated to produce an ensemble of possible precipitation scenarios for each rainstorm event.The model consists of three key components, described further below: (a) time-varying precipitation distributions; (b) spatiotemporal noise generation; and (c) precipitation ensemble simulation.

Time-Varying Precipitation Distributions
Parametric precipitation distributions conditioned on large-scale atmospheric variables were modeled at each 0.03°grid cell.The distributions consist of two components.
1.A precipitation occurrence model that predicts the probability of precipitation occurrence P wet in a grid cell at time t using logistic regression (Sperandei, 2014): where PR m and PW m are ERA5/CESM2 precipitation and precipitable water linearly interpolated to the grid, (i.e., the conditioned large-scale atmospheric variables).For model fitting and validation, PR m and PW m were obtained from ERA5 rainstorms (see Section 3.3.1).For model simulation, PR m and PW m were obtained from CESM2 rainstorms, under the assumption that the statistical relationships derived from ERA5 and AORC hold for CESM2 (see Section 3.3.2). α 0 , α 1 , and α 2 are fitted parameters.2. A precipitation magnitude model that predicts the distribution of precipitation amount in a "wet" grid cell using a transformed nonstationary gamma distribution (TNGD).Let x represent the precipitation amount at a single grid cell.We assumed x follows a generalized gamma distribution with shape parameters a > 0, c > 0, and scale b > 0 (Stacy, 1962): Applying a transformation y = x c gives y ∼ Gamma (a',b') , where a′ = a and b′ = b c (see Text S3 in Supporting Information S1 for derivations): The gamma distribution has mean μ = a′b′ and variance σ 2 = a′ b′) 2 .To incorporate nonstationarity, the mean μ was modeled as a function of large-scale atmospheric variables (Scheuerer & Hamill, 2015): where μ c is the mean of a stationary gamma distribution fitted to all the ERA5 rainstorms' AORC precipitation data at the grid, and β 0 ,β 1 ,β 2 and β 3 are fitted parameters.The variance σ 2 and shape parameter c were kept constant over time.The parameter c in the TNGD controls the tail heaviness of the distribution, allowing better fitting of extreme precipitation values compared to a basic gamma distribution.The transformation to variable Y also creates a simple relationship between the mean and distribution parameters that simplifies parameter estimation.

Spatiotemporal Noise Generation
We followed the STREAM (Hartke et al., 2022) approach, which itself is based heavily on Pysteps (Nerini et al., 2017;Pulkkinen et al., 2019), to simulate spatiotemporally correlated Gaussian noise that replicates the spatial and temporal dependency structure of the actual precipitation.Given an AORC precipitation field at the first time step t 0 of a rainstorm, we divided the field into overlapping spatial windows and computed the 2D Fourier amplitude spectrum within each window (Figure 2a).The amplitude spectra capture the spatial correlation of precipitation within each window while omitting phase information that determines the exact location of precipitation.For windows with less than 10% precipitation coverage, the global amplitude spectrum of the entire field was used.
Next, we filtered Gaussian white noise locally using the amplitude spectrum from each window: where n p 1 ,p 2 is the spatially correlated noise field, |R p 1 ,p 2 | is the local amplitude spectrum for the window at the position (p 1 ,p 2 ), N is the Fourier Transform of uncorrelated Gaussian white noise, and (∘) is pointwise multiplication.This generates noise with local anisotropy and spatial dependence matching that of the original AORC precipitation within the window.The final nonstationary noise field is obtained by summing the local noise from all windows.We chose an empirical window size of (128, 128) grid cells to balance capturing localized precipitation structures and maintaining sufficient samples within each window.A 0.3 overlap ratio between windows was chosen to improve the smoothness of the final noise field.This means that when delineating the windows along the horizontal direction, 30% of each window overlapped with the previous window.The same overlapping strategy was also applied in the vertical direction.For any overlapping regions between windows, the noise values were averaged.
The correlated noise field at the next time step t + 1 was modeled by a lag-1 autoregressive (AR(1)) model (Figure 2b): where n(t + 1,i,j) is the noise value to be calculated at time t + 1 and grid cell position (i,j), n(t,i v t ,j u t ) is the noise value at time t, which means the noise is advected by the 850 hPa north-south and east-west wind vector (v t , u t ) from position (i v t ,j u t ) at time t to position (i,j) at time t + 1, ñ(t + 1,i,j) is the new spatially correlated noise generated at t + 1, and ρ is the correlation coefficient between AORC precipitation at t and t + 1.The wind vectors were obtained from ERA5/CESM2 850 hPa wind fields.
Generating space-time correlated noise requires amplitude spectra from high-resolution AORC precipitation.For CESM2 rainstorms, corresponding high-resolution precipitation fields are not available.To obtain plausible spatial precipitation structures, we used a precipitation field matching method to find surrogate high-resolution precipitation fields from historical AORC data (Figure 3).This matching assumes rainstorms with similar large-scale atmospheric conditions will exhibit similar precipitation patterns and spatial correlations.We defined the distance between the CESM2 atmospheric variable fields at time step t and the ERA5 atmospheric variable fields at another time step t' as: where N is the number of grid cells along the horizontal direction in the CESM2/ERA5 atmospheric variable fields, M is the number of grid cells along the vertical direction, PR ERA5,i,j and PR CESM2,i,j are the precipitation values at grid cell position (i, j) from ERA5 and CESM2 fields, respectively.These atmospheric variables were min-max normalized before calculation.For each CESM2 time step t, we computed the distance d between it and every time step t' from the ERA5 rainstorms during 1979-2021.We then selected the five ERA5 time steps with the smallest distance to t.Finally, we randomly sampled one of those five time steps, with a probability proportional to the inverse of its distance to t.This ensures that ERA5 time steps with more similar atmospheric variables to the CESM2 time step have a higher probability of being sampled.The matching is performed independently at each time step.To account for seasonality, only the ERA5 time steps from the same or adjacent months to the current CESM2 time step were considered in calculating the distance.For example, if the CESM2 time step is in March, the pool of ERA5 time steps is constrained to February through April.The sampled time step provides an AORC precipitation field serving as the surrogate high-resolution precipitation, whose amplitude spectra are then used to generate spatially correlated noise.In effect, this finds historical precipitation fields that resemble the spatial structure of a given CESM2 rainstorm.It only uses local and global amplitude spectra and discards the phase information, thus reducing effects from location mismatch between AORC and CESM2 precipitation.There can be some rare cases where a local window in the domain contained CESM2 precipitation but lacked precipitation in the corresponding region of the sampled AORC field.In such cases, there is not enough high-resolution AORC precipitation data in that window to compute a local amplitude spectrum.As an approximation, we used the global amplitude spectrum, computed from the entire sampled AORC field, as the local amplitude for that window.Since the sampled AORC fields may lack temporal continuity, the correlation coefficient ρ in Equation 6was computed based on CESM2 precipitation when simulating CESM2 rainstorms.

Precipitation Ensemble Simulation
The actual precipitation field was obtained by transforming the noise field through the time-varying parametric precipitation distributions described in Section 3.2.1 (Figure 2c).The correlated noise field was first converted to a probability field using the CDF of the standard Gaussian distribution.If the transformed probability P(t) at one grid is lower than the current dry probability P dry (t) = 1 P wet (t), the simulated precipitation was set to zero.If P (t) exceeded P dry (t), the excess was rescaled to 0-1 as P′(t) = P(t) P dry (t) 1 P dry (t) .This P′(t) was then converted to y(t) using the inverse CDF of the fitted nonstationary Gamma distribution (Equation 3).The final precipitation x(t) was calculated as x(t) = y(t) 1/c .Multiple noise realizations for a rainstorm can be transformed to generate an ensemble of precipitation realizations.Each realization represents a possible precipitation scenario consistent with the current large-scale atmospheric state.

Experimental Design
An overview of the model experiments is as follows: The SRG was first fitted based on AORC precipitation and ERA5 large-scale atmospheric variables and then cross-validated (Section 3.3.1).Then, the SRG was used to simulate precipitation fields based on large-scale atmospheric variables from CESM2 rainstorms (Section 3.3.2).This assumes that the ERA5 and bias-corrected CESM2 rainstorms share similar characteristics, such that the established ERA5-AORC relationships can be applied to downscale CESM2 storms.The model experiments were conducted using computational resources provided by the Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison (Center for High Throughput Computing, 2006).

Model Fitting and Evaluation
The precipitation occurrence and magnitude models were fitted for each month from December-May using ERA5 rainstorms for 1979-2021.Fitting was performed at each 0.03°grid in the study domain (29-50°N, 79-113°W), consisting of 1,024 × 630 grids.6-hr AORC precipitation data and associated large-scale atmospheric variables (ERA5 precipitation and precipitable water) were extracted from the identified rainstorms.For logistic regression, AORC data were converted to 0 (dry) or 1 (wet) using a 0.2 mm threshold.The Python Scikit-learn package (Pedregosa et al., 2011) was used to estimate the logistic regression parameters using the maximum likelihood method.For TNGD, samples with AORC precipitation >0.2 mm were first fitted to a stationary generalized gamma model using Scipy's method of moments (Virtanen et al., 2020).The estimated shape parameter c transformed the precipitation data as y = x c .The new variable y was then fitted to the nonstationary gamma distribution (Equation 3).The distribution parameters were estimated by minimizing the mean continuous ranked probability score (CRPS, units in mm, see Hersbach, 2000) using the optimization function in Scipy: where F t is the CDF of the nonstationary gamma distribution at time t, y t is the transformed AORC precipitation, n is the sample size, and 1(y t ≤ z) is a Heaviside step function, which is 1 if y t ≤ z and 0 otherwise.CRPS measures the distance between CDFs, so minimizing the mean CRPS minimizes the difference between the parametric and empirical CDFs of observations at the grid.Fitting the model by each month was intended to account for seasonal effect.Although large-scale atmospheric variables can incorporate seasonal variations, statistical relationships between local precipitation and atmospheric variables can change across months or seasons.Estimating monthly varying TNGD parameters thereby incorporates such seasonal variations.However, monthly parameter estimation reduces the sample size for each month and may influence fitting quality.We would recommend users to consider both sample size and fitting performance when deciding whether to separate the fitting by month or season.
We used five-fold cross-validation to evaluate out-of-sample model performance.The original 43 years of data were split into five folds.Each fold used 80% of the data for fitting and 20% for simulation and validation.To simulate an ERA5-based rainstorm, we first generated correlated noise fields based on AORC precipitation.The noise fields were converted to precipitation through the precipitation distributions conditioned on ERA5 largescale atmospheric variables.We simulated 20 realizations per rainstorm, representing possible precipitation patterns other than the AORC precipitation.Looping through the five folds generated out-of-sample simulations for all ERA5 rainstorms.Key rainstorm characteristics, including total precipitation, rainstorm area, and spacetime autocorrelation, were compared between simulations and AORC data.The Brier Score (BS) was used to assess the accuracy of predicted precipitation occurrence (wet or dry) at each grid (Benedetti, 2010): where K is the number of realizations, x t, k is the kth simulated precipitation at time t, and x t is the AORC precipitation.BS ranges from 0 to 1, with 0 for perfect accuracy and 1 for highly inaccurate predictions.We also computed the Reduction CRPS as RCRPS( Ft ,x t ) = CRPS( Ft ,x t )/σ 0 to assess the accuracy of simulated precipitation distribution at each grid cell (Trinh et al., 2013), where Ft is the empirical CDF of the 20 precipitation realizations at time t and σ 0 is the climatological standard deviation of the AORC precipitation.Compared to the original CRPS, RCRPS is normalized and thus independent of the precipitation magnitude at each grid, allowing spatial comparison of model performance across grids and subbasins.A lower value of RCRPS indicates a smaller difference between the CDFs of simulated and AORC precipitation.

CESM2-Based Rainstorm Simulations
Random precipitation fields were generated for storms identified from 10 CESM2 ensemble members.For a given event, surrogate high-resolution precipitation fields were sampled using the matching method described above (Section 3.2.2).Correlated noise was then generated and transformed to precipitation based on the precipitation distributions conditioning on CESM2 large-scale atmospheric variables.To reduce computation, we simulated a single realization for each rainstorm in the 10 100-year members.The resulting storm total precipitation distributions were then compared against those from the AORC data.To estimate extreme precipitation frequency, we selected simulated rainstorms with annual maximum 1-day, 3-day, 7-day, and storm total precipitation and generated 20 additional realizations for each event.Bootstrap resampling was implemented to estimate return levels of annual maxima: (a) In each bootstrap sample, we randomly resampled 1,000 years with replacement and selected one of the 20 annual maxima realizations for each year.(b) The empirical CDF of annual maximum precipitation was estimated from this bootstrap sample using the Weibull plotting position.(c) Steps 1 and 2 were repeated 5,000 times to estimate mean return levels and their uncertainties.

ERA5-Based Rainstorm Identification and Simulation
The STARCH method identified intense IVT events and associated rainstorms over the basin (e.g., Figures 1 and  4a-4c).The tracked IVT objects exhibited an elongated narrow structure, with IVT intensities gradually decreasing from the center toward the edge.In most cases, only one IVT object was present at a time, moving steadily from southwest to northeast across the basin.Previous studies have linked strong water vapor transport with heavy precipitation in the basin (Lavers & Villarini, 2013, 2015;Nakamura et al., 2013;Steinschneider & Lall, 2016).These IVT events were frequently associated with extratropical cyclones or fronts (e.g., Figure 1c, Nayak & Villarini, 2017;Ralph et al., 2018), but they can also be caused by other mechanisms such as Great Plain Low-level Jet, Rossby waves, or tropical cyclones (Cook et al., 2008;Gori et al., 2022;Moore et al., 2019).Good agreement was found between ERA5 precipitation patterns and AORC data, suggesting that ERA5 simulates realistic large-scale atmospheric environments in which the rainstorms occurred.However, precipitation locations and spatial distributions still differed at local scales, suggesting that local precipitation is partly determined by mesoscale or finer-scale processes that are not fully resolved by ERA5.
Figure 4 shows a rainstorm in December 2019, with typical atmospheric conditions conducive to heavy precipitation: an extratropical cyclone system with intense IVT and high precipitable water covering the basin (Figures 4a-4c).The Gaussian noise field captured local spatial structures of the AORC precipitation (Figure 4d).For example, the noise in the northwest of the domain was smooth and continuous in the first time step, while the lower-right was fragmented, consistent with the AORC precipitation pattern (Figure 4g).The simulated spatiotemporal precipitation patterns were similar to the AORC data, exhibiting a comma-like structure (Figures 4e-4g).The model simulated heavy precipitation in the regions of high ERA5 precipitation and precipitable water, showing positive relationships between the mean of precipitation distributions and large-scale atmospheric variables.Notably, however, the local patterns and hotspots of simulated precipitation differed from the AORC and varied across realizations, reflecting the stochasticity in small-scale precipitation processes that could occur under the same large-scale atmospheric conditions.The noise fields in Figure 4d exhibited strong irregularity and local variations, which appeared to be largely reduced in the final precipitation simulations.This is because a large portion of the irregular noise outside the rainstorm was converted to zero precipitation based on the precipitation occurrence probability.In contrast, noise values within the rainstorm area directly inherited the corresponding spatial structures from the AORC precipitation fields, which resulted in smoother spatial variations and more regular patterns.
The simulated and AORC total precipitation fields for three extreme rainstorms preceding the 1997, 2011, and 2018 Mississippi Floods are shown in Figure 5.The rainstorms exhibited similar spatial orientations and coverage, with intense precipitation concentrated over the Lower Mississippi and Ohio-Tennessee subbasins.The simulated precipitation patterns and amounts generally agreed with the AORC data.For instance, total precipitation realizations for the 2018 event ranged from 47 to 57 mm with a mean of 51 mm, close to the 50 mm amount from AORC.Intense precipitation locations and shapes, however, differed from the AORC precipitation since the simulations were intended to represent multiple possible extreme scenarios.The AORC precipitation displayed stronger anisotropy and long-range spatial dependence (e.g., the more elongated arc of intense precipitation in the 1997 event, Figure 5a).Also, the AORC precipitation exhibited more continuous spatial variations, while the simulations contained more fragmented, noise-like textures.This suggests some discrepancies remained between the simulated and AORC precipitation structures that were not fully captured by the rainfall generator.
Consistency also existed in the temporal evolution of total precipitation, with the SRG ensemble spread (i.e., min to max) of 20 realizations containing most of the AORC time series (Figure 6).Notably, the SRG ensemble mean only represents the general trend of simulations and does not necessarily match the AORC data.This is because the AORC precipitation is also a single realization of an underlying random process, and therefore may not always locate at the mean at every time point.The spatial autocorrelation functions (ACF) of simulated rainstorms were generally consistent with the AORC data, though slightly lower from 0 to 400 km (Figures 7a and 7b).The rainstorms exhibited short-term dependence over time, with temporal ACF dropping below 0.2 after the 12-hr lag (Figure 7c), supporting the AR(1) temporal noise model.The lag-1 ACF coefficients matched the AORC, except being lower for the 2018 event.The spatiotemporal patterns, precipitation time series, and ACFs were also examined for the remaining ERA5 rainstorms, yielding similar agreement.
Basin-average total precipitation and rainstorm areas were compared between the SRG ensemble mean and AORC for each rainstorm.Despite the model generating variant spatiotemporal patterns across realizations, the average total precipitation was highly correlated with the AORC data (correlation ρ = 0.97, Figure 8a).This consistency is partly because the time-varying precipitation distributions were conditioned on ERA5 large-scale atmospheric variables, which indirectly limits the range of possible precipitation values.Another explanation stems from the atmospheric water balance: the net water vapor flux entering the river basin-governed by largescale atmospheric variables-should balance the total precipitation during an event, thus constraining possible precipitation amounts.There was also agreement in total and intense (>20 mm) precipitation areas between simulations and reference (ρ = 0.90 and 0.96, Figures 8b and 8c).For smaller rainstorms under 1.5 million km 2 , however, the model tended to overestimate total precipitation area.This overestimation may result from the coarse ERA5 data not fully resolving the AORC precipitation patterns and overpredicting precipitation  occurrence in some cases.Figure 1b at 18 hr, for example, illustrated an example where ERA5 generated greater precipitation coverage than the AORC data.
Model performance was also evaluated at each grid cell using the Brier Score and RCRPS.The average BS over the basin was 0.11, suggesting good accuracy in simulating precipitation occurrence (Figure 9a).Comparatively higher BS was seen in the Ohio-Tennessee (an average of 0.14), Lower Mississippi (0.13), and Upper Mississippi (0.12) subbasins (Figure 9b).The average RCRPS was 0.18 over the basin, indicating satisfactory agreement between simulated and AORC precipitation distributions (Figure 9c).The Ohio-Tennessee subbasin had the  highest average RCRPS of 0.25, followed by 0.22 in the Lower Mississippi and 0.2 in the Upper Mississippi subbasins (Figure 9d).The BS and RCRPS patterns suggest higher SRG uncertainties in regions experiencing more extreme and frequent rainstorms.A potential explanation is that these regions exhibit more complex precipitation processes, for example, localized convection or orographic lift, leading to more uncertain relationships between large-scale atmospheric environments and local precipitation (see Barlow et al., 2019;Schumacher & Johnson, 2006;J. A. Smith & Baeck, 2015;Su et al., 2023).

CESM2-Based Rainstorm Identification and Simulation
Similar rainstorm characteristics were found between the rainstorms identified from the 1979-2021 ERA5 data and those from the bias-corrected CESM2 in the same period (Figure 10).The average annual number of rainstorms from December-May was both 25 in ERA5 and CESM2 (Figure 10a).Average storm duration was around 50 hr for both data sets, with maximum durations reaching 240 hr in ERA5 and 258 hr in CESM2 (Figure 10b).The average precipitation area over the basin was about 1.3 million km 2 , approximately 40% of the basin area (Figure 10c).The average precipitation rate of the CESM2 rainstorms was 1% lower than that of the ERA5 rainstorms, while average precipitable water and IVT were 3% higher (Figures 10d-10f).However, CESM2 rainstorms had a larger standard deviation in average precipitation rate (16% higher) and lower standard deviations in average precipitable water and IVT (8.6% and 8.0% lower).CESM2 rainstorms also exhibited more extreme values, with the largest average precipitation rate of 8.8 mm/6-hr compared to 7 mm/6-hr in ERA5 rainstorms.Overall, bias-corrected CESM2 data produced rainstorms generally consistent with those from the ERA5 data, allowing them to be used for precipitation simulations.However, more complex biases may still exist in the CESM2 data after correction which can introduce some uncertainty into the simulated results.
Spatiotemporal patterns of a CESM2 rainstorm in January 1973 are shown in Figure 11.The matched AORC precipitation, despite being sampled from different time steps, exhibited similar locations and shapes to the CESM2 precipitation, providing reasonable approximations of precipitation spatial structures.The simulated precipitation fields more closely resembled the CESM2 patterns rather than the sampled AORC fields.This is because precipitation occurrence and magnitude were mainly determined by precipitation distributions conditioned on CESM2 large-scale atmospheric variables, while only the spatial correlation was drawn from the AORC data.The SRG realizations had very similar precipitation coverage and shape, with differences primarily in intense precipitation locations.This indicates the SRG ensemble members were not fully independent due to their shared reliance on the same large-scale atmospheric conditions.Two additional simulated CESM2 rainstorms are also provided in Figures S3-S6 of Supporting Information S1.
One realization of simulated total precipitation and associated atmospheric variables for three extreme rainstorms are shown in Figure 12.The events produced basin-average precipitation of 67, 75, and 91 mm, roughly equivalent to 50-year, 100-year, and 1,000-year events based on the GEV distribution fitted to the reference data (GEV fitting shown in Figure 15a).All three rainstorms occurred in early spring and were characterized by longlived (>160 hr) water vapor transport from the Gulf of Mexico.Consistent with ERA5 rainstorms, the simulated precipitation patterns exhibited southwest-northeast orientations, with intense precipitation near high CESM2 precipitation and precipitable water.The 1,000-year event displayed the greatest area of heavy precipitation (>300 mm), covering nearly the entire Lower Mississippi and parts of adjacent subbasins.Multiple peaks occurred in the simulated precipitation series (Figure 13), which is also an evident trait in ERA5 rainstorms (Figure 6).Despite showing similar trends to the original CESM2 precipitation, SRG ensemble members showed greater variability in temporal precipitation series.This enhanced variability in the SRG ensembles is more evident in the total precipitation amounts.For example, the total precipitation produced by the SRG ranged from 79 to 99 mm for the April 1962 event, whereas the bias-corrected CESM2 data showed 78 mm (Figure 13c).This suggests that the SRG introduces variability not only in the spatiotemporal precipitation patterns but also in the total precipitation amount of a rainstorm event.This randomness could potentially cover more extreme scenarios and lead to important variability in flood responses.
Since the simulated CESM2 rainstorms incorporated climate change signals, we separated the SRG simulations into historical  and future (2022-2050) periods to compare with AORC data.The simulated total precipitation distributions of CESM2 rainstorms from historical period (10,941 events) generally agreed with the AORC data for each month (Figure 14).An exception was May, where the model tended to simulate more events with 5-10 mm total precipitation.The simulated distribution of total precipitation from December-May also agreed with the AORC data (see Figure S7 in Supporting Information S1).The L-kurtosis, a measure of tail heaviness based on L-moments (T.J. Smith et al., 2023), was computed for each month using 1,000 bootstraps (Table 1).The average L-kurtosis was 0.153 for CESM2-based simulations during 1979-2021 and 0.167 for the  AORC, suggesting simulated rainstorms exhibited similar but slightly lighter tail behavior.However, the simulations exhibited significantly lighter tails (lower L-kurtosis) in December and heavier tails (higher L-kurtosis) in May, likely resulting from remaining biases in CESM2 data.Noticeable shifts in the precipitation distributions can be seen between future and historical periods of CESM2-based simulations (Figure 14).This suggests a greater mean and wider spread in total precipitation distribution in future climate, with increased probability density in the extremes.
Figure 15a compares the storm total annual precipitation maxima from CESM2-based rainstorms in the historical period to those from the AORC data.The SRG simulations showed good agreement with the AORC annual maxima, especially in the tail of the distribution.The simulated annual maxima were slightly lower than AORC data at around 10-year return period, possibly due to sampling uncertainty in the limited AORC records.Figures 15b-15d compare the annual maximum 1-day, 3-day, and 7-day precipitation between CESM2-based rainstorms and AORC data during 1979-2021, also showing overall agreement and consistent tail behavior.
The SRG tends to simulate larger 1-day and 3-day precipitation than AORC between 1-to 10-year return periods, which can be attributed to the remaining bias in the CESM2 data that are difficult to remove.In contrast, despite the original bias-corrected CESM2 precipitation producing 1-day and 3-day annual maxima close to the AORC data, it significantly underestimated extreme 7-day and storm total precipitation.This suggests that CESM2 data still contained significant bias in temporal and spatial correlation structures and may not be suitable for directly representing extreme scenarios, especially for long-duration events.These biases have been substantially reduced by the SRG, as shown in Figures 15a and 15d.
For comparison, GEV distributions were fitted to the AORC annual maxima in Figure 15a-15d using L-moments (T.J. Smith et al., 2023).The GEV showed a closer fit to weak-to-moderate extreme events (1-10 years return periods) and increased faster than CESM2-based simulations after 100 years.This is especially evident for 1-day annual maximum precipitation, where the SRG simulations estimated a lower 500-year return level of 20 mm compared to 24 mm from the GEV distribution.This likely results from the GEV exhibiting a heavy, unbounded tail (shape parameter <0 except for 3-day precipitation).However, the GEV tail behavior can be highly uncertain due to fitting difficulty given the limited sample size.The tail heaviness of extreme precipitation has been investigated and discussed by many studies (e.g., Koutsoyiannis, 2004;Marani & Ignaccolo, 2015;National Research Council, 1994;Papalexiou & Koutsoyiannis, 2013;Serinaldi & Kilsby, 2014).It is widely recognized that commonly-used extreme value distributions (including GEV) may not adequately represent the underlying precipitation processes, causing high uncertainties in the tails.
Figure 15 also shows return level curves from SRG simulations in the future period ).An overall increase from 1-to 500-year precipitation can be seen for 1-day and 3-day annual maxima.For annual maximum 7-day and storm total precipitation, increase mostly occurred between 1-year and 100-year events and gradually reduced after 100 years.This possibly implies changes in temporal precipitation structures in long-duration extreme events-precipitation is more concentrated during the peak period of rainstorm, while the total depth has fewer changes.Since 10 ensemble members only provided a limited number of annual maxima for the historical and future periods, the tail of return level plots was heavily controlled by the largest one or two events and can possess significant sampling uncertainty.To better quantify future changes in very rare extreme events, it is recommended to include simulations from more ensemble members.
Using all the simulated annual maximum events from 1951 to 2050, we extended the return level curve to 1,000 years (Figure 16).The 1,000-year return levels for 1-day, 3-day, 7-day, and storm total precipitation were 22, 46, 74, and 86 mm, respectively.The extended return level curves exhibited seemingly bounded tails, as the increment in precipitation depth slows down with increasing return periods.This implies a maximum possible event in the basin, which is physically reasonable because the precipitation was constrained by large-scale atmospheric variables and space-time correlation structures.Figure 16 demonstrates the ability of StormLab to produce extremely rare events via hundreds or even thousands of years of simulations.However, including the entire period of 1951-2050 neglected climate differences between different periods and can underestimate extremes in the future climate.It is therefore recommended to select a 30-year or smaller window of simulations, similar to the separation used in Figure 15, to construct the return level curve.

Contribution to Design Floods in the Lower Mississippi River Basin
The simulations from our SRG potentially contribute to the design flood in the MRB by extending the choice of historical rainstorms to over 1,000 synthetic years of events based on CESM2 data.The design flood study of V. A. Myers (1959) revealed that it is a critical sequence of extreme rainstorms-as opposed to a single precipitation event-that is likely to cause the most extreme flood in the lower MRB (Su et al., 2023).In other words, the flood magnitude depends on storm arrival orders in the MRB and time of concentration.Flood discharge due to  downstream precipitation can be amplified by peak discharges coming from upstream rainstorms.While developing new hypothetical storm sequences is beyond the scope of this study, we highlight that rainstorm sequences can be determined directly based on their arrivals in CESM2 simulations, providing a large number of more physically reasonable sequences compared to "manual creation" of one (as in Hypo-Flood 58A) or several sequences.Note that a rainstorm sequence with 1,000-year total precipitation does not necessarily result in a flood with 1,000-year peak discharge.This is because flood responses are determined by many other factors, including watershed antecedent conditions, spatiotemporal precipitation patterns, and arrival locations of rainstorms.Therefore, future work will be to integrate the precipitation simulations with large-scale hydrologic models to analyze peak discharge outputs and identify rainstorm sequences that result in floods at certain average recurrence intervals.

Model Uncertainty
Multiple sources of uncertainties exist in the model that can influence the simulation results.This section discusses four major sources of uncertainty:

Uncertainty From Data Sources
Biases in data sets contributed to model uncertainties in different ways.The AORC data, created by downscaling daily NLDAS-2 and Stage IV precipitation, may not fully capture the extreme precipitation values at hourly scales, introducing bias when fitting precipitation distributions and extracting precipitation spatial correlation structures.Meanwhile, although ERA5 is corrected by extensive observations via data assimilation, biases in large-scale atmospheric variables can be higher in earlier periods (e.g., before 2,000) due to fewer available observational data.A key source of uncertainty arises from the way in which ERA5 and CESM2 data are used-while the rainfall generator was fitted based on AORC precipitation and ERA5 large-scale atmospheric variables, CESM2 data were used as covariates to simulate rainstorms, assuming the fitted statistical relationships can be transferred to CESM2.Though bias correction of CESM2 was used to reduce its biases relative to ERA5, this could not eliminate discrepancies in spatial and temporal correlation structures.For instance, CESM2 tends to generate smoother precipitation patterns and larger areas of intense IVT and precipitable water than ERA5.To reduce this uncertainty, more advanced bias correction methods could be applied to correct temporal and spatial structures of CESM2 data (e.g., Cannon, 2018;Mathieu;Robin & Vrac, 2021;Vrac & Thao, 2020).However, one must be careful, as such multi-aspect bias corrections may damage or lose information in the original CESM2 data, such as climate change signals.

Uncertainty From Rainstorm Identification
It can be challenging to accurately identify the start and end of individual rainstorm events.In principle, identifying and tracking IVT objects is easier than tracking precipitation patterns due to their relatively regular shape and movement.This is evidenced by the satisfactory tracking results shown in Figures 1, 4, and 11, where IVT events and associated rainstorms were correctly identified.Nevertheless, the STARCH parameters are often empirically determined based on the visual quality of tracking patterns, which can introduce uncertainty.For example, lowering the high IVT threshold from 500 to 250 kg m 1 s 1 allows for identifying more weak and moderate water vapor transport events that were excluded in the study.Also, since the IVT intensity is projected to increase under climate change, using constant threshold parameters may result in identifying more persistent and longer duration IVT events in the future time periods (see Zhao, 2020).Changing the STARCH parameters, especially the IVT thresholds, can impact the identified rainstorm samples and further affect distribution fitting and SRG simulations.Users need to adjust the parameter setting based on the specific application purposes, for example, to focus solely on extremes or to include all the precipitation events for long-term, continuous simulations.The main purposes of applying rainstorm identification were to isolate extreme events from the continuous precipitation record and to enhance SRG computation efficiency.This is because with the help of parallel computing, the SRG can simulate hundreds of individual rainstorms at the same time, significantly reducing the computation time.However, it should be noted that StormLab also supports long-term, continuous precipitation simulation without any rainstorm pre-selection or identification, albeit at higher computational cost.
The selection of tracking domain also depends on the purpose of SRG simulations.Since most of the water vapor is transported from the Gulf of Mexico, the current tracking domain covering the MRB is sufficient to identify relevant IVT events.If the SRG are used to simulate local precipitation events over a much smaller watershed, the tracking domain should be reduced accordingly.In this case, the tracking may also need to be implemented using finer-resolution precipitation and IVT data sets rather than the coarse CESM2 data.

Uncertainty From Precipitation Distributions and Covariate Selection
The logistic regression model achieved good prediction accuracy in precipitation occurrence (see Brier Scores in Figure 9), although spatial variations existed across sub-basins.Greater uncertainties originated from fitting the TNGD due to its more complicated structure and larger number of parameters (i.e., Equations 2 and 4).However, by examining the quantile-quantile plots, we found that the TNGD model achieved a good fit to the AORC precipitation over the study domain (see Text S4 in Supporting Information S1).This is likely because TNGD has improved fitting at distribution tails with the transformation parameter c, in comparison to fitting a simple gamma distribution.However, some grids still exhibited under-or overestimation of extreme precipitation (see Figures S9-S14 in Supporting Information S1).Also, the TNGD parameters fields contained evident rough and noisy textures (see Figures S1 and S2 in Supporting Information S1), likely because the TNGD was fitted individually for each grid cell without imposing any smoothness constraints across the domain.The irregular textures may also stem from the NLDAS and StageIV data sets used to create the AORC precipitation fields.In principle, the distribution parameters should vary more smoothly over the domain since most terrain does not change rapidly in the MRB.Potential solutions would be to add a spatial smoothing regularization term to the objective function (Equation 8) or to apply a smoothing kernel to the fitted TNGD parameters.
In the TNGD model, the transformed variable y was assumed to have a distribution mean being a non-linear function of large-scale atmospheric variables and a constant standard deviation (i.e., Equation 4).Using a constant variance and transformation parameter c aimed to simplify model fitting and structure, which has been applied in many previous studies (see Bechler et al., 2015;Gutmann et al., 2022;Newman et al., 2015).Given the limited length of data, estimating a nonstationary transformation parameter c (e.g., the shape parameter of the generalized gamma distribution) can introduce significant uncertainty.However, the actual variance of precipitation can increase with large-scale precipitation and precipitable water.Assuming a constant variance may underestimate precipitation variability under extreme atmospheric variables.One convenient solution is to scale the precipitation variance σ 2 (t) with its mean μ(t) (Hartke et al., 2022;Scheuerer & Hamill, 2015): where μ c and σ 2 c are the mean and variance of a stationary gamma distribution, and β 4 , β 5 are fitted scaling parameters.When β 5 = 1, the variance σ 2 (t) changes linearly with μ(t).StormLab can accommodate such forms of nonstationary variance or more complex relationships.More realistic and accurate statistical connections between local precipitation and large-scale atmospheric variables can significantly improve model performance, which is a common challenge in stochastic simulation and statistical downscaling.
Covariate selection is also an important source of uncertainty since extreme precipitation in the MRB is influenced by various factors with significant spatial heterogeneity.On the local scale, extreme precipitation at one grid cell is correlated with large-scale precipitation, precipitable water, and water vapor transport.Local precipitation can also be affected by vertical wind velocity and Convective Available Potential Energy (CAPE) (Kim et al., 2022;Z. Wang & Moyer, 2023).On regional or global scales, teleconnections and large-scale atmospheric circulations play an important role.Rossby wave breaking and blocking patterns have been linked to large-scale precipitation and floods over the central United States (Hirschboeck, 1987;Moore et al., 2019).Precipitation in different sub-basins also exhibits varying associations with climate patterns like El Niño-Southern Oscillation, Arctic Oscillation (AO), North Atlantic Oscillation, and Pacific North American Oscillation (PNA) (see Mallakpour & Villarini, 2016;Su et al., 2023;S.-Y. Wang et al., 2014;Wise et al., 2018).In essence, the relationships between local precipitation and large-scale atmospheric dynamics are complex and location-dependent.In the TNGD model, we simplified this relationship by only selecting large-scale precipitation and precipitable water due to their robust correlation with local precipitation.Also, the teleconnection impacts can be implicitly reflected in the large-scale atmospheric variables used in the TNGD model.However, eliminating other relevant factors like vertical velocity or CAPE can introduce uncertainty in representing precipitation variability.Therefore, one must carefully balance the information gained from each additional covariate and the risk of overfitting and collinearity.In this study, P and PW were chosen as a trade-off between model complexity and performance, resulting in satisfactory precipitation simulation over the basin.Improvements could be made by identifying key covariates in each sub-basin based on physical reasoning and their correlation strength with local precipitation.
Another key function of large-scale atmospheric covariates is to convey climate change signals.The TNGD formulation allows the distributions of local precipitation to evolve with P and PW in the future periods.However, while maximum precipitable water scales at approximately 7%/K following the Clausius-Clapeyron equation (Held & Soden, 2006), extreme precipitation generally does not follow the same scaling, depending on factors such as precipitation type, duration, and location (Fowler et al., 2021;Lochbihler et al., 2017;X. Zhang et al., 2017).This suggests that precipitable water, when used as the covariate, may increase faster than actual precipitation in future projections, resulting in overestimation of the distribution mean.Therefore, to accurately account for climate change impacts, one must carefully select covariates that exhibit similar scaling with temperature as extreme precipitation.A possible modification is to use an alternative covariate suggested by O'Gorman and Schneider (2009) that incorporates vertical velocity and moist-adiabatic lapse rate, which scales better with extreme precipitation under temperature changes.

Uncertainty From Precipitation Spatiotemporal Structures
While the SRG generated overall consistent rainstorm patterns and shapes as compared to the AORC precipitation, approximating precipitation structures using spatiotemporal Gaussian noise can introduce uncertainties in precipitation simulations.Although the local-window FFT improved representation of local precipitation spatial structures, it reduced samples in each window and decreased the accuracy of the amplitude spectrum that is related to spatial correlation.Precipitation fields in the local windows were assumed to be stationary, which can neglect some highly localized, anisotropic precipitation processes like local convection and orographic effects.Cropping the precipitation field by window also lost precipitation long-range structural information, weakening long-range dependence in simulated precipitation (see Section 4.1).Non-linear transformation via precipitation distributions can reduce correlation strength between Gaussian noise (Papalexiou, 2022), which may explain the reduced spatial ACF of simulated rainstorms (see Section 4.1 and Figure 7).One possible improvement is applying a residual spectrum to the original amplitude spectrum as compensation for correlation loss during transformation.For the time dimension, the AR(1) model omitted precipitation temporal dependence beyond lag 2 and can also be influenced by the quality and resolution of 850 hPa wind data, introducing additional uncertainties.
Another key source of uncertainty came from the precipitation field matching method for CESM2 rainstorms.The sampled AORC precipitation-although sharing similar large-scale atmospheric variable patterns-can have different spatial coverage and local precipitation structures to the unknown precipitation fields of CESM2 rainstorms.Also, independent sampling at each time step reduced temporal coherence of spatial structures, compared to ERA5 rainstorm simulations based on consecutive AORC fields.Since AORC data only covered 43 years, there might be no suitable AORC precipitation pattern for some CESM2 rainstorms, for example, rare events with a mixture of two or more trajectories of moisture transport in the basin.These uncertainties arose from the rainfall generator's reliance on high-resolution precipitation for spatial correlation information.To address this limitation, one potential improvement is to extend the high-resolution precipitation samples with simulations from convective-permitting regional climate models (e.g., Rasmussen et al., 2023).This allows the model to draw from a more diverse set of high-resolution precipitation patterns to match the CESM2 rainstorms.Additionally, one may use parametric spatial correlation models for each local window, with coefficients estimated based on current large-scale atmospheric variables.Another way is leveraging deep learning models like generative adversarial networks (e.g., pix2pix, Isola et al., 2017) or diffusion models (e.g., SR3, Saharia et al., 2022) to predict high-resolution precipitation or spectrum fields.

Summary and Conclusions
In this study, a SRG StormLab was proposed to simulate 6-hr, 0.03°rainstorm fields driven by water vapor transport in the MRB.The model used a local-window FFT-based approach to simulate nonstationary Gaussian noise fields that capture local spatial structures and temporal evolution of precipitation.The noise fields were transformed into actual precipitation amounts through precipitation distributions at each grid cell, where precipitation occurrence was modeled with logistic regression, and precipitation magnitude was represented by a transformed nonstationary Gamma distribution (TNGD).The precipitation probability and mean of TNGD were modeled as functions of large-scale atmospheric variables (i.e., ERA5/CESM2 precipitation and precipitable water) to incorporate spatial and temporal nonstationarity.Given the large-scale atmospheric variable fields of a rainstorm, the model could simulate multiple precipitation realizations representing possible precipitation scenarios.
Spatiotemporal rainstorms driven by strong IVT were identified from 1979 to 2021 based on ERA5 reanalysis and high-resolution AORC precipitation using a rainstorm tracking method STARCH.Precipitation distributions were fitted monthly at each grid cell for December-May based on the identified rainstorms.

Figure 2 .
Figure 2. Process of spatiotemporal precipitation generation.(a) Generating spatially correlated noise using local window fast Fourier transform based on high-resolution Analysis of Record for Calibration precipitation; (b) Generating temporally correlated noise using the AR(1) model.ρ Represents the correlation coefficient between precipitation fields at t and t + 1; (c) Transforming noise fields to precipitation fields through CDFs.

Figure 3 .
Figure 3. Process of precipitation field matching method.(a) Randomly selecting an ERA5 time step with the closest atmospheric variable fields to those at the target CESM2 time step; (b) Using Analysis of Record for Calibration precipitation at the sampled ERA5 time step to provide local (high-resolution) amplitude spectra for the precipitation simulation at the target CESM2 time step.

Figure 4 .
Figure 4. Spatiotemporal patterns of an ERA5-based rainstorm starting at 00:00 UTC on 29 December 2019.(a) ERA5 IVT; (b) ERA5 precipitation; (c) ERA5 precipitable water; (d) The first realization of noise; (e) The first realization of simulated precipitation; (f) The second realization of simulated precipitation.The two realizations are drawn from out-of-sample simulations; (g) Reference Analysis of Record for Calibration precipitation.For patterns of precipitation occurrence probability and transformed nonstationary gamma distribution parameters, please see Figure S1 in Supporting Information S1.

Figure 5 .
Figure 5.Total precipitation patterns from out-of-sample simulations (single realization) and Analysis of Record for Calibration for rainstorms in (a-b) February 1997, (c-d) April 2011, and (e-f) February 2018.

Figure 6 .
Figure 6.Time series of basin-average total precipitation from out-of-sample simulations (blue) and Analysis of Record for Calibration (gray) for rainstorms in (a) February 1997, (b) April 2011, and (c) February 2018.The stochastic rainfall generator ensemble contains 20 realizations.

Figure 8 .
Figure 8.Comparison of rainstorm characteristics between stochastic rainfall generator ensemble mean and Analysis of Record for Calibration for all ERA5 rainstorms.(a) Basin-average total precipitation; (b) Total precipitation area; (c) Total intense precipitation area (>20 mm).

Figure 9 .
Figure 9. (a) Brier Scores across the simulation domain; (b) Boxplots of Brier Scores in the Mississippi River Basin (MRB) and the five subbasins; (c) Reduction CRPSs across the simulation domain; (d) Boxplots of Reduction CRPSs in the MRB and the five subbasins.The five subbasins are also shown in (a), including (1) Lower Mississippi, (2) Arkansas-Red, (3) Missouri, (4) Upper Mississippi, and (5) Ohio-Tennessee subbasins.

Figure 10 .
Figure 10.Distributions of rainstorm characteristics from ERA5 (blue) and CESM2 (red) during 1979-2021.(a) Annual storm number from Dec-May; (b) Storm duration; (c) Average precipitation area; (d) Average precipitation rate; (e) Average precipitable water; (f) Average integrated water vapor flux.N is the total number of identified rainstorms.It should be noted that these characteristics were calculated directly from ERA5 and bias-corrected CESM2 data, without any simulation results from StormLab.

Figure 13 .
Figure 13.Time series of basin-average total precipitation of CESM2 rainstorms in (a) March 2022, (b) April 2029, and (c) April 1962.The dark blue line represents the rainstorm realization shown in Figure 12, and the gray lines represent the 20 stochastic rainfall generator ensemble members.The light blue line represents the biascorrected CESM2 data.

Figure 14 .
Figure 14.Distributions of basin-average total precipitation based on Analysis of Record for Calibration data (black) and simulated CESM2 rainstorms in the historical (1979-2021, solid blue line) and future (2022-2050, dashed blue line) periods.
(a) uncertainty from data sources; (b) uncertainty from rainstorm identification; (c) uncertainty from precipitation distributions and covariate selection; and (d) uncertainty from precipitation spatiotemporal structures.

Figure 15 .
Figure 15.Estimated return levels of annual maximum precipitation based on stochastic rainfall generator (SRG) simulations in the historical (1979-2021, dark blue line) and future (2022-2050, light blue line) periods.The red line represents the GEV distribution fitted to Analysis of Record for Calibration (AORC) annual maximum records.The dashed orange line represents the original bias-corrected CESM2 precipitation.The SRG simulations, CESM2, and AORC annual maxima were plotted based on empirical return periods estimated by Weibull plotting position.The uncertainty bound (shaded area) represents the 95% uncertainty interval (2.5th-97.5thpercentile) of the SRG simulations at each return period across the 5,000 bootstrap samples.(a) Storm total annual precipitation maxima; (b) 1-day annual precipitation maxima; (c) 3-day annual precipitation maxima; (d) 7-day annual precipitation maxima.

Figure 16 .
Figure 16.Estimated return levels of annual maximum precipitation based on stochastic rainfall generator (SRG) simulations during 1951-2050 from 10 CESM2 ensemble members.The empirical return periods were estimated based on Weibull plotting position.The uncertainty bound (shaded area) represents the 95% uncertainty interval (2.5th-97.5thpercentile) of the SRG simulations at each return period across the 5,000 bootstrap samples.
To evaluate the model LIU ET AL.