Advancing Space‐Time Simulation of Random Fields: From Storms to Cyclones and Beyond

Realistic stochastic simulation of hydro‐environmental fluxes in space and time, such as rainfall, is challenging yet of paramount importance to inform environmental risk analysis and decision making under uncertainty. Here, we advance random fields simulation by introducing the concepts of general velocity fields and general anisotropy transformations. This expands the capabilities of the so‐called Complete Stochastic Modeling Solution (CoSMoS) framework enabling the simulation of random fields (RF's) preserving: (a) any non‐Gaussian marginal distribution, (b) any spatiotemporal correlation structure (STCS), (c) general advection expressed by velocity fields with locally varying speed and direction, and (d) locally varying anisotropy. We also introduce new copula‐based STCS's and provide conditions guaranteeing their positive definiteness. To illustrate the potential of CoSMoS, we simulate RF's with complex patterns and motion mimicking rainfall storms moving across an area, spiraling fields resembling weather cyclones, fields converging to (or diverging from) a point, and colliding air masses. The proposed methodology is implemented in the freely available CoSMoS R package.

The four above-mentioned properties (i.e., skewness, dependence, anisotropy, and advection) are related and influence one another in several ways at different spatiotemporal scales. Synoptic weather systems, such as (extra) tropical cyclones, comprise large scale air masses rotating around a strong center of low atmospheric pressure exhibiting spiral-like spatial structure directly related to rotational motion. Much of the significant weather observed in association with such systems tends to be concentrated within narrow bands called frontal zones containing large mesoscale precipitation areas (rain bands) roughly oriented parallel to the front lines (Houze et al., 1976). Mesoscale precipitation areas are composed by regions of cumulus convective precipitation, commonly known as convective cells (Amorocho & Wu, 1977;Gupta & Waymire, 1979;Houze, 2014;Houze et al., 1976;McMurdie & Houze, 2006). Bands are nearly perpendicular to the direction of movement of the front, and the motion of these bands is sometimes faster than the motion of the general storm system (Amorocho & Wu, 1977;Ippolito, 1989). Cells within a rain band tend to move along or slightly ahead of the direction of the front, and the more intense cells tend to elongate in their direction of motion (Houze, 2018;Houze et al., 1990;Ippolito, 1989). This relationship between the shape of intense cells and advection direction was also recognized by Moszkowicz (2000), studying the spatial correlation functions of radar data for scales from two to tens of kilometers.
The interconnection among correlation, advection, and anisotropy plays a key role also in the quantification of these properties. Rainfall fields show an apparent motion resulting from the combined effect of the winds at some steering level and the systematic precipitation growth and dissipation (Germann et al., 2006). This suggests the possibility of analyzing the rainfall process on either Eulerian or Lagrangian coordinates, whereby the former refer to a fixed reference system (e.g., the ground), while the latter to a reference system moving with the precipitation field. Taking the motion into account, one can define Eulerian and Lagrangian variants of temporal correlation of radar rainfall fields (Zawadzki, 1973). Thus, the maximum of the cross-correlation function or the maximum of the lag-correlation structure can be used as recognizable features of the rainfall pattern. These features are then followed as they move in space across the study area (i.e., the Eulerian reference system), thus quantifying storm velocity and direction (Bacchi & Kottegoda, 1995;Niemczynowicz, 1987;Rinehart & Garvey, 1978;Shaw, 1983;Vischel et al., 2011;Zawadzki, 1973).
Concerning the link between anisotropy and correlation, anisotropy of geophysical fields is often studied by using two-dimensional spatial correlation functions, which allow for the analysis of both the scale over which patterns occur and the direction of the pattern (Niemi et al., 2014). Two-dimensional spatial correlation functions and the corresponding Fourier power spectra (where the latter are the Fourier transform of former according to the Wiener-Khinchin theorem) have been extensively used to study the spatial structure of precipitation fields (e.g., Cassiraga et al., 2020;Gyasi-Agyei, 2016;Krajewski, 1987;Mandapaka & Qin, 2013;Niemi et al., 2014;Sinclair & Pegram, 2005;Zawadzki, 1973). They are also the basis of methods, such as the Generalized Scale Invariance (Lovejoy & Schertzer, 1985), that allow one to quantify the scaling of anisotropic systems, thus accounting for the different anisotropy of cells and rain bands (see also Ramanathan & Satyanarayana, 2019).
The foregoing discussion highlights the importance of developing effective stochastic models for environmental/geophysical flows, such as rainfall, that account for possibly complex (non-Gaussian) marginal distributions, spatiotemporal correlation, anisotropy, and advection. However, building models incorporating these properties is challenging. Recent attempts concerning wind speed and lightning are provided for example by Gneiting et al. (2006), Youngman and Stephenson (2016), and North et al. (2020), while the problem of modeling rainfall was tackled, for example, by Paschalis et al. (2013), Leblois and Creutin (2013), Niemi et al. (2016), Nerini et al. (2017), as discussed below in more depth.
However, the existing methods generally focus on a subset of those characteristics, thus lacking some accuracy/generality in the description of the remaining features. Building on Papalexiou (2018) and Papalexiou and Serinaldi (2020), this study proposes a general framework that allows consistent modeling of all the foregoing properties along with extensions and additional features. We extend the Complete Stochastic Modeling Solution (CoSMoS) framework, which allows for an accurate reproduction of marginal distributions and spatiotemporal correlations, to incorporate general advection velocity fields and anisotropy. CoS-MoS is a general stochastic modeling framework that can be applied to any geophysical process, yet rainfall is the most natural example for this type of models.

Random Fields and Space-Time Correlations
Random fields (RF's) offer an integrated approach to model the spatial variability and temporal evolution of environmental processes. Natural processes are continuous in space and time, yet in practice RF's are typically implemented in discrete space and time. Let Ω be the sample space of a random experiment; a spatiotemporal RF is a stochastic process  In a plain language, in Lagrangian perspective, the observer follows an individual fluid parcel as it moves through space and time. Conversely, from Eulerian standpoint, the same observer focuses on specific locations in the space through which the fluid flows as time passes. In the context of rainfall fields, records of rain gauges are Eulerian point images of the ground-level rain fields produced by the moving bands (Amorocho & Wu, 1977), while radar echoes enable Lagrangian analysis by tracking the field motion in a sequence of radar images at subsequent time steps. Therefore, taking the motion into account, the Eulerian version of spatiotemporal correlation of radar fields depends on the motion of the storm and changes in structure that occur within it, while the Lagrangian variant is defined relative to the storm coordinates and therefore is independent of the storm motion, which is indeed filtered out in Equation 5 (Zawadzki, 1973). An important consequence of the equivalence in Equation 6 is that a Lagrangian STCS   does not peak at   0 for   0. A first description of this characteristic of Lagrangian covariance functions, which is called dimple effect, is due to Kent et al. (2011), while Cuevas et al. (2017 provided a discussion for the Gneiting class of STCFs. In essence, a STCS has a dimple if the rv  ( , ) s here now t is more correlated with  ( , ) s there then t than with  ( , ) s there now t .
Here we offer further insights on the dimple effect and show, for a constant advection vector, how the STCS varies depending on the motion direction in space. The STCS values in a Lagrangian framework depend not only on the distance and time lag between two rv's, but also on the direction of the field's movement. In the Cartesian plane the direction from point A to B is the angle of the line connecting the two points (see Figure A1). Thus, an STCS expressing an advected field depends on three independent variables (distance, time, and direction) and can be visualized in a 3D plot if a variable is fixed.
Let a space-time process, for example a storm, be observed over an  n n spatial grid comprising 2 n rv's, with advection velocity v x y x y v v , ( , ) ( , )   4 4 . Thus, the field is moving over the  45 diagonal covering a distance of 4 2 spatial units every time step. If the untransformed STCS  is the non-separable AMHW (Figure 1a), it can be modified based on Equation 6 to become a Lagrangian STCS   to account for any advection speed. In this case, the Lagrangian STCS shows a dimple in space and time (Figure 2 and Figure 3).
For temporal lags   0, dimple in space is apparent as the STCS peaks at nonzero distance ( Figure 2) in contrast to the static STCS ( Figure 1). If the advection vector is constant, the Lagrangian STCS peaks in the opposite direction of movement (Figures 2a and 2b), that is at  225 , for lag   1 and distance 4 2 (the distance the field covers in one time step). For temporal lag   2, it peaks at distance 8 2 (Figures 2c and 2d) since this is the distance the RF covers in two time steps. These peaks are the "global" maxima of the STCS PAPALEXIOU ET AL.  coincide (e.g.,   0To simulate Lagrangian Gaussian RF's we use multivariate (vector) autoregressive models of order p (MAR( ) p ) modifying the approach for static RF's simulation described in Papalexiou and Serinaldi (2020). Specifically: 1. An RF is described by a set of rv's , i j  over an  n n grid (see the example in Figure A2a in Appendix A). (1) D  in Figure A2d). are therefore plugged into the multivariate Yule-Walker equations (e.g., Lütkepohl, 2005) to estimate the MAR( ) p parameters. 5. The MAR( ) p model allows the simulation of vectors of 2 n rv's that preserve the desired STCS and the velocity vector.
quantile function of the rv  and Φ ( )   is the standard Gaussian cdf. Yet this transformation is non-linear and when applied to a bivariate Gaussian rv ( , ) due to the maximal property of the bivariate Gaussian distribution (Gebelein, 1941;Lancaster, 1957;Maung, 1941). Analytical expressions linking  , i j   and  , i j   exist only for specific cases including the uniform (Baum, 1957) and the Lognormal (Matalas, 1967) distributions. In the past, numerical or asymptotic approximations were used based on Hermite-Chebyshev polynomial expansions (Lancaster, 1957;van der Geest, 1998).
Here, we use the scheme introduced in Papalexiou (2018) and extended to STCS in Papalexiou and Serinaldi (2020). The link between the correlations of processes with Gaussian and non-Gaussian marginals is formed by using a parsimonious nonlinear function that acts as a correlation transformation function (CTF), that is, where  0 b and  0 c are parameters (estimation details are given in Papalexiou (2018), and software implementation in CoSMoS R package ). This function yields the val- Papalexiou and Serinaldi (2020) proposed two ways to apply the CTF; here, we use the second one, which assumes a well-defined STCS    ( , )  describing the Gaussian RF's, and defines the STCS of the target . This approach has two advantages: (a) it guarantees that the Gaussian RF's are described by positive definite STCS that enables their simulation; and (b), in the case of Lagrangian RF's discussed in the next section, it allows the characterization of the target process by the STCS of the parent non-Lagrangian Gaussian RF's. This overcomes the inconvenience of the STCS's of the target process, which are spatially varying, and are not described by simple STCS expressions in a Lagrangian coordinate system. PAPALEXIOU ET AL.  Therefore, hereafter, we characterize a space-time process (in its full extent) by four components: (a) the STCS of a parent Gaussian RF, which is a positive definite parametric STCS, (b) the marginal probability distribution of the target RF's, (c) the velocity vector field describing the advection, and (d) the anisotropy of the target RF's.
10.1029/2020WR029466 8 of 26 : Figure 5f demonstrates how a Gaussian value is transfromed to a non-Gaussian one. The arrows show a Gaussian z  0 8 . transformed to x  4 8 . . As anticipated, the values of the Gaussian RF's follow the standard Gaussian distribution (Figure 5e), while the transformed ones follow the desired BrXII distribution (red line in Figure 5g). 6. Figure 5h: the non-Gaussian RF's, transformed from the Gaussian ones (Figure 5d), have the desired STCS ( Figure 5b) and the desired marginal distribution (Figure 5g).

Lagrangian and Non-Gaussian Intermittent Random Fields
The previous non-Gaussian RF's ( Figure 5h) evolve in time without a systematic movement toward a specific direction, or else, their velocity vector is  , (0,0) x y v . We can generate Lagrangian non-Gaussian RF's by using Lagrangian Gaussian RF's with the aforementioned CoSMoS framework. In particular, the parent Gaussian RF's are simulated (see section 3) based on Lagrangian STCS's to account for advection velocity; plugging these RF's into the CoSMoS algorithm leads to RF's that move in time with constant velocity , Here we demonstrate simulation of intermittent Lagrangian RF's with mixed-type marginal distributions characterized by probability mass at zero 0 p and continuous distributions describing positive values. Fields of precipitation, wind and so forth, are typically intermittent in space and time, and thus, can be described by mixed-type distributions. The method of using any mixed-type marginal to generate time series with prescribed correlation structure was introduced in Papalexiou (2018) and generalized for RF's in Papalexiou and Serinaldi (2020). In particular, Papalexiou (2010) discussed the case of mixed-type marginals with probability mass at arbitrary points, while Papalexiou and Serinaldi (2020) highlighted the case of mixedtype marginals with probability mass at the boundaries of the sample space and continuous distribution in-between. Here, we recall the expressions of the mixed-type cdf where     0 denotes the conditional non-Gaussian rv which is defined in  (0, ) in this case, and u z  ( ) are probabilities corresponding to standard Gaussian values. Specifically, the mixed-type quantile function is used to transform the Gaussian RF's.

Advancing Random Fields Simulation
Over the next sections we advance spatiotemporal modeling of RF's by coupling them with velocity fields and general forms of anisotropy. We show that the advection velocity (defined as speed and direction of mass) can vary in space (or even in time) and can be described by a velocity field. A velocity field associates a vector describing the direction and the speed of the velocity to each point of space. We use these vectors to represent the advection velocity at each grid node of the discretized domain. Velocity fields are ubiquitous in nature and describe, for example, the movement of liquid or gas particles in space and time. Here, we focus on static velocity fields that could express for example a steady state system, yet we can also extend the simulations by using velocity fields changing in time. Our approach to couple RF with velocity fields is based on the idea that the velocity speed and direction at each point of the field modifies the STCS based on the Lagrangian distance   in Equation 6. In this case, and in contrast to constant velocity (fixed speed and direction), the Lagrangian distance   depends both on the direction between two points and their absolute coordinates; yet   is expressed by the same formula.
Next, we demonstrate RF's simulations coupled with velocity fields and anisotropy having the potential to model a wide-range of natural phenomena including liquid or gas motion. We show cases of RF's that rotate, have low and high speeds in different regions, converge from opposite directions, converge to a point or diverge from a point, and mimic cyclones, wavy patterns, and storms.  tropy. Our aim in the following demonstrations is to show the potential of the modeling framework and not to model actual cases based on observed data. Consistent and reliable estimation of advection, anisotropy, and STCS from real data is challenging and out of this study's scope. For brevity, we summarize the above four components for each experiment (Table 1) and avoid repeating details in the next sections.

Rotating Random Fields
The first demonstration regards rotating fields over a center point with velocity vectors described by uniform curl. The velocity vector at a point ( , ) x y , in an  n n field, is given by where x w and y w are scaling factors, and ( , ) ( / , / ) x y n n 0 0 2 2  is the center of the velocity field. The velocity speed for  x y w w is constant over fixed-radius circles; this leads to clockwise rotating fields with zero speed in the center and maximum speed in the corners (see velocity vectors in Figure 7a). Note, for  x y w w , the vectors' magnitudes are constant over fixed ellipses. Here, we simulate intermittent, rotating, non-Gaussian RF's (see Table 1 for details) described by the velocity field in Equation 10 for   1 / 3 x y w w . Inspection of a few frames (Figure 7) shows that the intermittent fields rotate clockwise; for example, the high intensity region spotted on the lower-left corner (Figure 7

Low-High Velocity Speed
A case of potential interest involves advection fields with low speed in a region and high speed in another. We demonstrate RF's with low-high speeds that are described by the velocity field    , , x y x y w y w x v For  x y w w , velocity vectors are symmetric with respect to the main diagonal  y x. The speed increases progressively from zero in the lower-left corner to its maximum in the upper-right corner (see vectors in Figure 8a). The velocity direction progressively changes from  90 in the lower-right corner to  45 in the upper-right corner, and from  0 in upper left to  45 in the upper right. We generate 100 RF's (see Table 1 for details) based on the previously described velocity field using w w x y    7 2 75 2 / (the highest advection speed is seven and spotted in the upper-right corner in an  75 75 grid). The simulation shows horizontal advection from left to right in the upper part (see e.g., the high intensity region in Figures 8.4-8.7), and vertical advection from the lower-right corner toward the top (see e.g., Figures 8.9-8.12). The almost zero speed in the lower-left corner and the high speed in the upper-right corner are clearly demonstrated in Movie S8.

Hyperbolic Velocity Fields
There are cases where mass (e.g., gas, liquid, etc.) moves from opposite directions toward a point (colliding) and must escape from it in different directions to preserve mass balance. Such a velocity field could be described by with ( , ) ( / , / ) x y n n 0 0 2 2  . For  x y w w the velocity vectors are symmetric over the diagonals  y x and   y x. This case could fit to mass moving from northwest and southeast corners toward the center; mass speeds cancel out in the center, and, to preserve mass balance (if accumulation is not allowed), mass leaves the region from the southwest and northeast corners (see velocity vectors in Figure 9a). The advection "streamlines" are hyperbolas and their conjugates; the maximum speed, observed in the corners, is v n n w n ,  2 for an  n n field and   x y w w w . We simulate fields (see Table 1) coupled with this parabolic velocity field (Figure 9a of the simulated fields ( Figure 9) reveals the advection over the parabolas; yet the movement is clearly shown in Movie S9.

Radial Velocity Fields: Divergence and Convergence
Velocity fields with positive or negative divergence are omnipresent in nature and relate with movement of liquid or gas. For example, heated air expands toward all directions and the movement of air mass forms a velocity field. In such fields, there is a point that acts like a source as the flux moves outwards from it in PAPALEXIOU ET AL.

10.1029/2020WR029466
13 of 26  all directions, or else, the velocity vectors point outwards as shown in Figure 10a. In physics this is called positive divergence and the velocity vectors can be described by Here the divergence point is set in the center of the spatial domain, that is, ( , ) ( / , / ) x y n n 0 0 2 2  . For  x y w w the velocity vectors are symmetric over the center and have equal speed over circles with same radius (Figure 10a). Here we demonstrate generation of RF's having positive divergence for w w x y   6 2 75 / . A characteristic sequence of simulated RF's ( Figures 10.1-10.5) shows the divergent advection. A region of low intensity in the center of the field, surrounded by clusters of higher intensity (Figure 10.1), expands outwards over the next frames pushing the high intensity clusters out of the field (see Movie S10).
In contrast to the previous demonstration, velocity fields of gas or liquid can have negative divergence, or else, converge to a point. For example, cooling air contracts; the local decrease of the gas volume due to cooling combined with the external pressure, forces the gas particles to move inward or toward the low-pressure region (negative divergence). In general, such velocity fields can occur when a point acts like a sink causing flux to move toward it from all directions (Figure 11a). Velocity vectors with negative divergence can be described by PAPALEXIOU ET AL.

10.1029/2020WR029466
14 of 26 (1-11) a sequence of generated fields (Movie S10 clearly shows the expansion). Similarly to positive divergence, velocity vectors with negative divergence for  x y w w are symmetric over the convergence point and have constant speed over fixed-radius circles (Figure 11a). Here, we simulate convergent RF's for  x y w w having the same properties as the divergent fields (see Table 1). Note that in the Supporting Information we report general and compact representation of the velocity fields' equations with the examples presented here being special cases. A representative sequence of generated RF's  shows the negative divergent advection. A few clusters of high intensity, surrounding a large low-intensity region in the center of the field (Figure 11.45) move progressively toward the center decreasing eventually its extent (e.g., Figure 11.49). The convergence is better demonstrated in Movie S11.

Affine Anisotropy
Referring to Allard et al. (2016) for a detailed review of anisotropy in geostatistics, we can identify three elementary models called geometric, zonal, and separable, which can be combined to provide more complex anisotropies (Journel & Froidevaux, 1982). Among these elementary models, the geometric one is probably the most used in literature as discussed in section 1. It involves an affine transformation of the Cartesian coordinate space according to the relationship where  x and  y scale the coordinates of a point  ( , ) x y s T , and  is the angle of counterclockwise rotation around the origin (e.g., Chilès & Delfiner, 2009).
In the previous section, we demonstrated RF's with advection described by velocity fields. Locally varying advection can be combined with anisotropy. Here, we combine affine anisotropy with the velocity field described in Equation 13 with ( , ) ( , ) x y shows that the  75 75 generated RF's follow the prescribed velocity field. We clarify that Figure 12b shows how a regular rectangular grid is transformed after applying affine anisotropy (Equation 15) to its coordinates (see also Figures A2c and A2d for an intuitive explanation).

Simulation Mimicking Cyclones
The rigid "stretching-and-rotation" geometric transformation can be generalized in the spirit of velocity fields by introducing locally varying scaling and rotation factors, that is, using deformation fields based on local affine transformations (Ligas et al., 2019). This method enables the simulation of fields with convenient anisotropy mimicking for instance cyclonic shapes via a swirl-like coordinate transformation given by x y from the center of the deformation ( , ) x y 0 0 ,  is a rotation angle, and  is a scaling parameter controlling the swirl extension. Figure 13b provides an ex- where 1 w and 2 w are weight factors for each one of the two velocity fields. This velocity field (Figure 13a), when combined with swirl anisotropy (Figure 13b), could mimic the behavior of cyclones.
Here, we demonstrate a non-Gaussian (see Table 1) cyclone-mimicking simulation. The velocity field (Figure 13a) is specified by w 1 2 2 15 5  /( ) and w 2 2 15 5  /( ) so the maximum speed in the corner is five, and the swirl anisotropy has parameters    3 2 / and  18 b (Figure 13b). The simulation shows clear cyclonic advection and swirl anisotropy (Figure 13). The RF's rotate but at the same time masses move toward the center as it is clearly shown in Movie S13.

Wavy Anisotropy
We demonstrate the potential of coupling advection with new forms of anisotropy. We experiment here by creating a form of wavy anisotropy described by where the parameters  1 ,  2 ,  3 , and  control the characteristics of the "wave." Wavy patterns are omnipresent in nature and are observed in sand and soil, wood texture, rock formations, surface of liquids, etc. Figure 14b provides an example of this wavy deformation by showing a transformed rectangular grid of ( , ) x y points to ( , )   x y points using   0,   1 2,   2 1 / 2,    Paschalis et al. (2013) proposed a rainfall generator called STREAP (Space-Time Realizations of Areal Precipitation), which is based on the so-called "String of Beads" model (Pegram & Clothier, 2001), and forms one of the components of the so-called Advanced WEather GENerator for a two-dimensional grid (AWE-GEN-2 days; Peleg et al., 2017). STREAP assumes lognormal marginal distributions and relies on the simulation of latent Gaussian RF's to account for spatiotemporal correlation. As the RF's are simulated by Fast Fourier Transform (FFT), FFT symmetries are exploited to fold the simulated fields and mimic advection. The model seems not to incorporate anisotropy, FFT generation requires the simulation of much larger fields to avoid numerical artifacts due to field symmetries, and marginal back-transformation from Gaussian to Lognormal space introduces bias in the STCS's. Leblois and Creutin (2013) proposed a rainfall generator called SAMPO (Simulation of Advected Mesoscale Precipitations and their Occurrence), which relies on the simulation of latent Gaussian RF's similar to STREAP. In this case, latent fields are generated by a three-dimensional turning band method, which allows the introduction of spatiotemporal correlation, space-time anisotropy (to fulfill the so-called Taylor hypothesis (Taylor, 1938)), and advection (by using the concept of back-trajectories). Unlike STREAP, SAMPO reproduces the target spatio-temporal correlation by inflating the correlation structure of the latent Gaussian process via an anamorphosis function based on standardized Chebyshev-Hermite polynomials. SAMPO yields rainfall fields with inverse-Gaussian marginal distribution, which is skewed but asymptotically exponential, whereas rainfall tends to exhibit subexponential (heavy) tails at daily and subdaily scales of interest (Cavanaugh et al., 2015;Moustakis et al., 2021;Nerantzaki & Papalexiou, 2019;Papalexiou et al., 2013;Papalexiou & Koutsoyiannis, 2016;Rajulapati et al., 2020;Serinaldi & Kilsby, 2014). Similar to STREAP, SAMPO requires the simulation of fields covering larger areas to simulate advection effects. Benoit et al. (2018) proposed a model conceptually similar to SAMPO, in which the field generation is based on Cholesky decomposition and sequential simulation, while the skewed distribution of positive rainfall is approximated by a power transformation of latent Gaussian fields.

Existing Space-Time Rainfall Models Featuring Advection and Anisotropy
Latent Gaussian fields are also the core of the model proposed by Niemi et al. (2016), which generalized the so-called STEPS model (Short Term Ensemble Prediction System; Bowler et al., 2004;Seed et al., 2013) introducing scale-varying anisotropy by GSI (hereinafter STEPS-GSI). Similar to STREAP, the STEPS-GSI model assumes lognormal distribution for positive rainfall, and it requires the generation of larger fields (as it uses FFT simulation) as well as a post-processing of simulated rainfall to adjust field mean and standard deviation. The model does not explicitly address the problem of correlation bias due to logarithmic back-transformation of the latent approximately Gaussian fields. This problem also affects the nonstationary stochastic rainfall generator based on the short-space Fourier transform (SSFT) described by Nerini et al. (2017). The SSFT-based model enables performing a local moving window Fourier analysis to account for changes in the anisotropy and extent of correlation structures across a field. This approach differs from STEPS-GSI, which accounts for changes across scales (e.g., local cell structures within rain bands) instead of changes across a geographic area. Similar to STREAP and STEPS-GSI, the SSFT-based model assumes approximately lognormal distribution for positive rainfall, requires post-processing of simulated fields to adjust field mean and standard deviation, and does not incorporate advection.

Storm Simulation
The concepts introduced in the previous sections can be used to advance stochastic simulation of storms and weather systems. Combining such concepts results in a framework that offers several advantages and advances: 1. It allows generating storm fields described by any non-Gaussian distribution and STCS preserving intermittency. 2. It is scale free, that is, it can reproduce the characteristics of storms at any temporal and spatial scale given that these characteristics are specified. 3. It allows general advection expressed by velocity fields with locally varying speed and direction. 4. It supports the use of general anisotropies beyond affine transformations.

10.1029/2020WR029466
18 of 26 5. It enables the generation of any number of RF's since the scheme is not based on pre-generated large fields limiting the length of simulation, and 6. since the proposed method relies on multivariate AR models, it is free from the constraint of generating RF's over a square region (  n n grids). Indeed, the procedure generates only rv's at specified spatial coordinates, which can define a regular or irregular grid as well as the nodes of a gauging network.
Here we demonstrate storm simulation by progressively adding more components in the stochastic model, showing three characteristic cases.
Constant-velocity, isotropic storms. Over large regions, storms might have varying velocity (both in speed and direction). Yet in small regions, constant velocity could be a valid assumption. Additionally, while radar observations typically show some form of affine anisotropy there are cases that storms can have isotropic patterns. Here, we mimic an isotropic and constant velocity storm by generating intermittent non-Gaussian RF's (see Table 1) with constant advection velocity  , (7,7) x y v . The strong space-time correlation is apparent in the simulation (Figure 15; Movie S15) and enables clusters of high intensity cells to persist in time, and intermittency (zero and non-zero regions) to persist in space. A characteristic cluster of high intensity cells ) reveals the north-eastward advection velocity entering the region from lower-left corner and exiting from the upper-right corner (see also Movie S15).
Constant-velocity, anisotropic storms. As mentioned above, most storm radar observations at fine spatiotemporal scales (e.g., 1 km and 10 min), or even satellite observations at larger spatial scales, indicate that storm cells or precipitation patterns exhibit some form of affine anisotropy (see section 7.1). Assuming constant velocity, we can easily add affine anisotropy in the simulated intermittent RF's. Here, we simulate a constant-velocity and anisotropic storm by adding affine anisotropy to the previous simulation with anisotropy parameters   2 x ,   1 y , and     /4 (see Equation 15). The anisotropy is apparent in the generated RF's ( Figure 16; Movie S16) and clearly observed in the cluster of high intensity cells. The selected values of  x and  y create elliptical anisotropy with the ratio of major to minor axis equal to two, while the rotation parameter  makes the major anisotropy axis perpendicular to the advection velocity, which is  , (7,7) x y v in this case.  coupled with velocity fields (with or without anisotropy) to simulate storms. Keeping the same STCS and marginal distribution as in previous two examples, we demonstrate a storm simulation with varying-velocity and anisotropy. In this case, the affine anisotropy parameters are set to  x  1 2 / ,   2 y , and     /4 placing the major axis over the  y x diagonal (see Figure 17b showing how a regular square grid is transformed). In contrast to the prior examples, where velocity fields had constant speed and direction, we create a velocity field with constant speed and varying direction (Figure 17a). In general, this is achieved by dividing the horizontal (   100 100 storm fields characterized by the Ali-Mikhail-Haq-Weibull STCS, BrXII distribution, movement from south-west to north east, and affine anisotropy (see Movie S16).

Varying
Here, we set y v x n , and   8 x y w w . This velocity field has constant speed at all points and curved "streamlines" with curvature and directions that vary over the region (Figure 17a). The generated RF's ( Figure 17; Movie S17) show the anticipated affine anisotropy and the movement of the fields over the "streamlines." The advection is more clearly visualized in Movie S17.

Discussion and Conclusions
This study advances the simulation of spatiotemporal RF's providing more realistic representations of rainfall, weather systems, or other hydro-environmental fluxes. We generalized the CoSMoS framework, which allows the simulation of RF's with specified non-Gaussian marginal distributions and STCS's. The proposed generalization is based on the concepts of general velocity fields and general anisotropy transformations. The former allows the introduction of locally varying advection, thus extending the commonly used uniform motion across an area or pure random motion. On the other hand, general anisotropy transformations yield space deformations going beyond the simple affine anisotropy transformation. Combining these two concepts enables the simulation of RF's with complex structures/patterns and motion across a specified spatial domain, mimicking for instance spiraling vortices, sink or source effects, or colliding air masses. As proof of concept, we provided a set of simulations introducing progressively additional elements of complexity. We also proposed a new copula-based spatiotemporal correlation functions enabling great flexibility in modeling spatiotemporal dependence, and provide theoretical conditions for their validity, that is, their positive definiteness.
As the saying goes "Rome wasn't built in a day"; this study focuses on model development and simulation strategies and offers a step forward toward a general and powerful stochastic modeling framework for environmental applications. Several topics require further investigation. Specifically, model inference, a crucial topic for real world applications, is not tackled here. As mentioned in the introduction, links among correlation, anisotropy, and advection, can be exploited to estimate the required parameters. Existing methods should be assessed in detail to form an estimation procedure that is consistent with the framework presented here; this explains why we did not report real-world examples.
Spatial variability of marginal distributions and STCS's (i.e., spatial non-stationarity) is another feature that can improve, for instance, the modeling of low-high pressure regions, and orographic effects on weather systems. Similarly, temporal non-stationarity that would allow deterministic evolution of the statistical properties of the RF's was not considered. Temporal evolution of velocity fields is another aspect related to real-world applications that can improve the simulation of storm patterns or wind fields, for example. Lastly, simulation at synoptic scales (more than 1,000 km) and very fine spatiotemporal resolution (e.g., 1 km) is challenging due to computational limitations. Some of these features are technically easy to implement, while others need extensive effort.
Extensions enabling the simulation over non-planar surfaces, such as quadratic surfaces (spheres and hyperboloids), are also relevant to open the door to stochastic simulation of geophysical fluxes at global scale. Furthermore, non-Gaussian dependence structures (e.g., using the t copula) can also be considered to account for quantile-varying spatial correlation, and upper tail dependence, thus better representing the different joint behavior of moderate and extreme values.
Finally, since this study extends the capabilities of CoSMoS, one of its natural applications is to generalize the DiPMaC algorithm (Precise Temporal Disaggregation Preserving Marginals and Correlations; ) that breaks down coarse-scale time series into finer scales. The improved field-compliant CoSMoS represents the core of an enhanced version of DiPMaC to allow for downscaling weather products available at coarse spatiotemporal resolutions, such as Earth System Model and Regional Climate Model outputs, or gridded reanalysis and satellite products. All previous aspects are under investigation and will be the topic of future communications.
where   0 and  denote scale and shape parameters. The BrXII has two shape parameters  