Multifractal properties of embedded convective structures in orographic precipitation : toward subgrid-scale predictability

Rain and cloud fields produced by fully nonlinear idealized cloud resolving numerical simulations of orographic convective precipitation display statistical multiscaling behavior, implying that multifractal diagnostics should provide a physically robust basis for the downscaling and sub-grid scale parameterizations of moist processes. Our results show that the horizontal scaling exponent function (and respective multiscaling parameters) of the simulated rainfall and cloud fields varies with atmospheric and terrain properties, particularly small-scale terrain spectra, atmospheric stability, and advective timescale. This implies that multifractal diagnostics of moist processes for these simulations are fundamentally transient, exhibiting complex nonlinear behavior depending on atmospheric conditions and terrain forcing at each location. A particularly robust behavior found here is the transition of the multifractal parameters between table and unstable cases, which has a clear physical correspondence to the transition from stratiform to organized (banded and cellular) convective regime. This result is reinforced by a similar behavior in the horizontal spectral exponent. Finally, our results indicate that although nonlinearly coupled fields (such as rain and clouds) have different scaling exponent functions, there are robust relationships with physical underpinnings between the scaling parameters that can be explored for hybrid dynamical-statistical downscaling.


Introduction
It is often observed that air masses impinging upon a mountain range lead to the development of an orographic cloud with shallow embedded convective structures, which change the rainfall pattern and amount considerably, and can lead to localized extreme values of rainfall.These localized extremes present a great challenge for forecasters and are responsible for mountain hazards including landslides, debris flows and flash floods.
The notion that embedded convective structures are the result of exponential growth of small-scale disturbances in an unstable stratified layer (Kuo, 1963) motivated the application of linear stability analyses to gain insight on the dynamics of orographic convective precipitation, including the estimation of the unstable growth rate, ω, for each small-scale disturbance mode of wavenumber k = (k x , k y , k z ) (Fuhrer and Schär, 2005): − c l m ln (T ) + g where N m is the moist Brunt-Väisälä frequency of the saturated layer (as defined by Emanuel, 1994), which has a negative value in the moist unstable case.In this inviscid formulation, the linear model predicts maximum growth for infinitesimally small horizontal wavelengths.Admitting that smallscale disturbances are ubiquitous in the real atmosphere, this does not agree with the observed finite spacing between bands observed in nature, e.g., in western Kyushu in Japan (Yoshizaki et al., 2000), in the Cévennes region in southern France (Miniscloux et al., 2001;Cosma et al., 2002) and the Oregon Coastal Range (Kirshbaum and Durran, 2005b;Kirshbaum et al., 2007b).Kirshbaum et al. (2007a) developed a Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
new (also inviscid) linear model that includes an upstream stable region where stationary small-scale lee waves are triggered and subsequently used as the initial disturbances for the convective region.Their model agrees with the results from numerical simulations that showed that such lee waves play a dominant role in triggering and organizing banded convection with finite spacing (Kirshbaum and Durran, 2005a;Fuhrer and Schär, 2007).Using both linear stability analysis and idealized numerical simulations, additional factors were also found to be important such as (i) cloud depth (which determines the minimum value of k z allowed in Eq. 1) (Kirshbaum and Durran, 2004); (ii) the advective timescale, i.e., the in-cloud residence time of air parcels which can be controlled by mean wind speed and mountain width (Furher and Schär, 2005); (iii) dry stability outside the cloud layer (Kirshbaum and Durran 2005a); and (iv) cloud base height (Kirshbaum et al., 2007a).
A different approach to the problem has gradually developed concurrent but separately in the last thirty years, where the (statistical) scaling behavior of physical processes is explored.In these processes the statistical properties of a field at different scales are related by a scale-changing operation (generally a power law) that involves only the scale ratio and a scaling exponent, the simplicity of which is very appealing for statistical downscaling applications.In geophysical fields it is usually found that the scaling is determined not by one, but by an infinity of scaling exponents given by the scaling exponent function, i.e., they are multiscaling.Schertzer and Lovejoy (1987) proposed a functional form for this function, the universal multifractal (UM) model, briefly presented in Sect. 4.There are a number of publications reporting multiscaling behavior in various geophysical fields, including cloud and rain fields, over various temporal (see e.g., de Lima and de Lima, 2009 for a review) and spatial ranges, e.g., Schertzer and Lovejoy (1987).Nykanen (2008) found scaling on horizontal maps of radar rain reflectivity; Gupta and Waymire (1993) found scaling on spatial fields of oceanic rainfall; Tessier et al. (1993) and Barros et al. (2004) found scaling in satellite-based fields of cloud radiances.Recently, Lovejoy et al. (2008a) looked at TRMM product 2A25 near surface reflectivity from rain measurements and found evidence of multiplicative cascades from planetary scales down to a few kilometers (∼ 4 km, the data set resolution) in the ensemble averaged scaling statistics.Also, note that the 2A25 product is an orbit based product with limited dimensions in the direction perpendicular to the satellite movement.They argued that the UM model with well posed fixed values of the scaling parameters is a good approach when considering the entire range of scales.Stolle et al. (2009) showed that data from global models (including ERA40 reanalysis and NOAA GFS) also accurately follow cascade statistics from nearly 20 000 km down to around 100 km, where there was a cut off by hyper viscosity of the numerical model linked to grid resolution.
However, there is no consensus on the specific values of the scaling exponents (and associated UM parameters) and scaling ranges in the literature, although some of the differences might be attributed to sample size or differences in measurements (sensors, measurement methods or measurement errors, etc.).Another view is that the scaling exponents vary depending on contextual environmental conditions, and therefore the observed variations can be attributed in part to physical processes.Over and Gupta (1994) also used 2-D oceanic rainfall fields from GATE (GARP Atlantic Tropical Experiment), and proposed a relation between a scaling parameter and large-scale average rain (argued to be an index of synoptic-scale weather conditions); Perica and Foufoula-Georgiou (1996) related scaling parameters computed from horizontal rainfall radar fields to the Convective Available Potential Energy (CAPE) of the pre-storm environment over the same area; Deidda (2000) and Deidda et al. (2004) also found a dependency in one of their scaling parameters computed from radar reflectivity on the large-scale mean rainfall rate; and more recently Nykanen (2008) linked the variation in the values of the multiscale statistical parameters computed from radar rainfall fields to the underlying topographic elevation, predominant orographic forcing and storm location, and movement relative to the orographic cross sections.These results suggest that UM parameters vary from event to event, or case study to case study as a function of atmospheric conditions, and local terrain parameters.If this behavior can be translated into robust relationships between UM parameters of particular states (e.g., precipitation) and the specific physical processes that determine the space-time evolution of these states (e.g., moist convection), then inference and induction modeling approaches can be used to develop quantitative dynamical models of the state of interest relying on knowledge of the specific physical processes alone.For instance, to predict precipitation fields over a wide range of spatial scales knowing atmospheric stability conditions and precipitation at one (coarser) scale.This type of modeling could be used to formulate parameterizations of unresolved processes in numerical weather and climate prediction models, and to predict states at spatial resolutions much finer than those resolved by models or observations, i.e., downscaling.In practical terms, this requires investigating possible relationships between physical properties and statistical parameters that would allow statistical downscaling from coarse to fine resolutions.
In the present work scaling analysis is performed on the output of highly idealized numerical simulations, in particular on surface accumulated rainfall fields and horizontal cross sections of cloud water mixing ratios at several vertical levels.This gives us the opportunity to investigate horizontal scaling at several vertical levels at the same time, obtaining a horizontal scaling regime for 3-D physically consistent fields, instead of 2-D fields used previously in the literature (see above).The behavior of the scaling parameters for varying configurations of small-scale terrain and upstream The scaling in different fields is compared in Sect.4.4.In Sect.4.5, scaling analysis based on Fourier power spectra is performed on the 3-D cloud fields.The final section (Sect.5) presents a summary and discussion of the main results.

Numerical simulations
The numerical model Advanced Weather and Research Forecasting (WRF-ARW) version 3.1 (Skamarock et al., 2008) is used to perform idealized cloud resolving simulations of conditionally unstable moist flow impinging upon a Gaussian shaped ridge (defined by Eq. 3), elongated on the cross-flow direction, y, and with a maximum terrain height h m = 600 m.When the flow is lifted by orography, it saturates and releases latent heat generating an unstable orographic cloud where the embedded convective structures will develop.Decay parameters σ x and σ y are 12.5 and 5 km, respectively, and the center of the ridge is defined at (x c , y c ) = (80, 60) km.The resulting modified Gaussian ridge is represented in Fig. 1a.In all simulations the horizontal resolution is 250 m and there are 70 vertical terrain-following levels, unequally distributed over the 15 km of the domain, with vertical spacing stretching from lower to higher levels in order to have better resolution in the region where orographic clouds will develop.A Rayleigh damping layer is introduced in the top 5 km to reduce spurious reflections at the top.To simplify the problem, the parameterizations for radiation or surface and boundary layer are turned off, and the effects of earth's rotation are also neglected.The 3-D Smagorinsky scheme is used for sub-grid turbulence closure.These highly idealized simulations are based on previous works of Kirshbaum and Durran (2005a, b) and Kirshbaum et al. (2007a, b) where the authors showed that they reproduce essential features of observed convective bands in the Oregon Coastal Range.The model resolution is fine enough to resolve the convective scale (see e.g., Fuhrer and Schär, 2007), and thus it is assumed that these simulations represent the appropriate convective structures.However, there are two important differences between the present simulations and previous work: (i) the Thompson microphysics scheme (available in WRF 3.1) including cold microphysics is used in this study instead of warm rain schemes, which may be difficult to justify in the case of orographic precipitation simulations; disturbances in the domain reinforcing certain wavelengths, and would not allow the upstream stable flow to flow around the ridge, eventually forcing it all to transpose the orographic barrier.
The upstream profiles are also highly idealized, initially being horizontally homogeneous.In order to constrain the parameter space, they are defined in a very simple manner by a background flow in the x direction (V = W = 0) without shear (U (z) = const), a constant stable dry Brunt-Väisälä frequency (N d = 0.01 s −1 ), a constant surface temperature, T s and a relative humidity profile defined by three separate layers: RH = 90 % from the surface to 2500 m, followed by linear decay from 90 % to 1 % from 2500 m to 3500 m, and then RH = 1 % above 3500 m.Small-scale topography is added by superposing a 2-D sinusoidal field with a maximum amplitude of 50 m (which is less than 10 % of h m ), and the same wavelength in the x and y directions (l sst x = l sst y = l sst ) (Fig. 1b).The base case simulation, called Sst10km, has l sst = 10 km, T s = 285 K and U = 10 ms −1 .A population of simulations is built from the base case by varying one of the simulation parameters (l sst , T s and U ) at a time, while keeping the other parameters constant (see Table 1 for a summary of the population of WRF simulations conducted in this study).The particular choice of parameters derives from the physical intuition given by linear stability analysis on embedded convection briefly presented above.In this simple setup N 2 m is varied by varying T s and the advective timescale, τ adv , is varied by varying U .The control simulation CTL has no small-scale terrain.

Dynamical interpretation: linear stability analysis
Simulation CTL shows weak and relatively disorganized banded structures with a spacing between them on the order of the minimum wavelength that can be represented by the model horizontal resolution (Fig. 2a), consistent with the inviscid linear model in the absence of small-scale terrain.
Here, the small-scale disturbances required for convective triggering should be introduced either by physical mechanisms (e.g., adjustment of the initial profiles to the largescale terrain), or nonphysical numerical errors.This is not a particularly relevant case for the atmosphere where smallscale roughness is always present at the surface, triggering lee waves that will play a dominant effect in embedded convection as discussed above.This fact becomes clear by comparison of CTL with simulations that have small-scale terrain, the latter showing intense and well-organized convective rain bands aligned parallel to the mean wind (Fig. 2b, c,  d), with significant precipitation enhancement: total domain precipitation can be more than doubled and maximum local rainfall intensity can increase more than 17 times, depending on l sst (Table 2).These results are qualitatively similar to Kirshbaum et al. (2007a) linear model predictions.However, for a single scale of terrain variability the analytical model predicts a band spacing equal to l sst whereas our simulations show banding with l sst /2 spacing.The reason is that while the linear model assumes that the upstream cloud edge position is fixed, the numerical simulations show that bands can form at different x positions, generated by terrain features that are 90 • out of phase, a nonlinear effect not captured by the linear model.
Although the idealized analytical models give interesting insights on governing parameters and physical mechanisms, they fail to predict the correct pattern except in particular idealized cases, and they are unable to produce realistic quantitative results.Besides different cloud edge positions, the linear models are also unable to capture other nonlinear effects observed in the numerical simulations such as decay of convective structures, different intensities between bands in the same simulation, and band narrowing (see Kirshbaum et al., 2007a).Furthermore, the linear models should only be valid (at most) in the initial stage of convective growth, while disturbance amplitudes are small.Thus, nonlinear models are necessary to improve our knowledge of embedded convection, and for quantitative predictions.In the next section, one particularly promising nonlinear model based on cascade scaling models from turbulence is introduced.

The universal multifractal model
The existence of scale invariance in atmospheric fields is often investigated by analyzing their statistical moments, which are expected to obey the generic multiscaling relation: where the angled brackets represent the statistical average, q is the moment order generalized to any positive real number, λ = L 0 l is the scale ratio, L 0 being the outer scale of the cascade (the largest scale of variability) and l the scale of the observation, or simulation.ϕ λ is a quantity that is on average conserved from scale to scale in a similar way to what happens in turbulence cascade models.The turbulent flux is non-dimensionalized and normalized, such that ϕ λ = 1.In general geophysical fields are multifractals, i.e. the exponent K is a nonlinear function of q, and thus an infinity of exponents are required to characterize the scaling behavior.However, one can use the UM framework to model the variation of K with q for a conserved process, reducing the problem to two parameters (α and C 1 ; Schertzer and Lovejoy, 1987): The Levy index, α, defined in the interval [0,2], indicates the degree of multifractality (α = 0 for monofractals); and the co-dimension of the mean singularity, C 1 , describes the sparseness or non-homogeneity of the mean of the process.
The moment order must be positive (q > 0) in this framework.There is no physical basis to assume that a general observable, f , (e.g., wind, temperature, rainfall, etc.) should be conserved on a scale by scale basis.In analogy with turbulence models, we may expect that observable fluctuations, f (l), over a distance, l, are related to a general turbulent conserved flux, ϕ l , by A third universal multifractal parameter is introduced in Eq. ( 6): the Hurst exponent (or nonconservation parameter), H .In a conservative process ϕ λ = cte, and thus H = 0 and K(1) = 0.But in the general case, the mean has a dependence on scale, < f >∼ l H .The other exponent, η, depends on the particular field being analyzed.Note that Eq. ( 6) is a generalization of the classical laws of turbulence.For example, in the Kolmogorov (1941) law for turbulent velocity fluctuations: v(l) = ε 1/3 l 1/3 , we have that H = 1/3, η = 1/3 and ϕ l = ε.In this case the observable (turbulent www.nonlin-processes-geophys.net/20/605/2013/ Nonlin.Processes Geophys., 20, 605-620, 2013 wind field, v) is not the direct result of a multiplicative cascade that role being reserved for the (cascade conserved) energy flux, ε.The cascade conserved flux can have more complex forms that include nonlinear interactions of different fluxes, e.g., the Corsin-Obukhov law for passive scalar advection (Obukhov, 1949;Corsin, 1951): ρ(l) = ξ 1/3 l 1/3 , where ξ = ϕ l = χ 3/2 ε −1/2 , χ being the scalar variance flux.
In these classic turbulence cases, the laws relating the observable fluctuations and conserved turbulent fluxes can be obtained from the governing dynamical equations, assuming isotropy, and via dimensional analysis.For rain and cloud processes (which generally cannot be assumed to behave as passive scalars), the values of η and H are unknown, and so is the physical nature of the conserved flux.However it is still possible to perform scaling analysis.We start by taking the simplifying assumption η = 1, as usually done in similar previous works, noticing that if ϕ λ is a pure multiplicative cascade then that is also true for ϕ η λ , although with different C 1 values (e.g., Tessier et al., 1993;Stolle et al., 2009).The next step is to remove the influence of the l H term.One method used for this purpose is to rely on the absolute value of a finite difference gradient to compute the fluctuations at the highest resolution available, f (l res ) (see e.g., Lavallée et al., 1993;Tessier et al., 1993): The flux estimates given by Eq. ( 7) are systematically degraded to lower resolutions (lower values of λ) by spatial averaging.The Double Trace Moment (DTM) technique (e.g., Lavallée et al., 1993) can be used subsequently to estimate α and C 1 .Notice that the values of the UM parameters estimated from data generated by numerical models are not the same as the ones obtained from observations due to cascade cutoff by model resolution and model domain, but it is possible to convert between them (Stolle et al., 2009).Nevertheless, incorrect representation of physical processes in the models cannot be corrected.In the present work we aim only to investigate the variation of the UM parameters in numerical simulations with similar setup, so no conversion is required.
Spectral analysis can also be used to investigate the behavior of a field over a wide range of scales.The spatial Fourier power spectrum, E(k x , k y ), is computed by multiplying the 2-D fast Fourier transform of a field by its complex conjugate, where k x and k y are the wavenumber components.The power spectrum is then averaged angularly about k x , k y = (0, 0) to yield what is usually called isotropic power spectrum E(k), with k = k 2 x + k 2 y .Spatial scaling invariance manifests itself as log-log linearity of the power spectrum in space: where β is the spectral exponent and the addition of −1 in the exponent is required due radial averaging is phase space (e.g., Turcotte, 1992;Lovejoy et al., 2008b).Least square regression is used to estimate β from a log-log plot of E(k) against k.
The spectral exponent is related to the UM parameters by (Tessier et al., 1993) Finally, we estimate H from the slope of the log-log plot of the first order structure function against the lag, as described by Harris et al. (2001), Nykanen ( 2008) and Lovejoy et al. (2008a, b).

Scaling in simulated 2-D fields
The empirical scaling analysis described above is applied to surface accumulated precipitation spatial fields over 5 h simulations, P 5h ac .Here the constant value 80 km (scaling analysis domain length) was used for the cascade outer scale, L 0 .Changing the chosen value of L 0 has no effect on the estimated value of the slope of the fitted lines, i.e., on K(q) and respective UM parameter estimation.Figure 3a-c show log-log plots of M q against λ for P 5h ac fields for three different simulations.The lines are only fitted to the range of scales that are well resolved by the model, i.e., larger than 5 x(= 1.25 km) and smaller than about L 0 /4(= 20 km).For q ≤ 2 the linear fits are good implying that statistical scale invariance predicted by Eq. ( 4) is a good representation, and that there are no characteristic scales or scale breaks for the considered ranges of q and λ (from ∼ 1 to ∼ 20 km).This range of q is similar to the ones found in literature for similar analysis (e.g Douglas and Barros, 2003;Lovejoy et al., 2008a, b;Nykanen, 2008;de Lima and de Lima, 2009).Some simulations show persistent deviations from linear behavior at particular scales, but the relation between these scales and l sst varies between simulations, suggesting that the smallscale terrain wavelength is not (statistically) a characteristic scale in the rain fields.The same spatial scaling analysis was repeated for the P 5h ac fields resulting from all simulations in Table 1, all of them showing the same robust scaling behavior (linear relations).
The scaling exponent function, K(q), and respective UM parameters can be used to quantify the multifractal behavior.Figure 4 shows that for q < 2 there is a very robust match between the empirically estimated scaling exponent function, K(q), (i.e., slopes of the fitted lines in log-log plots of M q against λ) from P 5h ac fields for the different simulations and the UM model curves (solid lines), computed from Eq. ( 5) with the UM parameters obtained from DTM technique, which is consistent with the literature (e.g., Nykanen, 2008;Lovejoy et al. 2008a;Stolle et al., 2009, among others).It is also important to notice the variation of K(q) with varying simulation attributes (l sst , T s and τ adv ).This variation is quantified in Fig. 5  Relationship between M q and λ for different values of q.M q is computed from P 5h ac from simulation (a) Sst5km, (b) Sst10km, (c) Sst15km.The scales considered on the analysis are marked by "o" while "x" marks the ones not considered.The plots are on a log-log scale.
which display a complex nonlinear variation with simulation attributes.A particularly remarkable feature is the abrupt transition in both UM parameters for T s between 275 K and 285 K.The transition has a physical correspondence when interpreted against the spatial rainfall patterns (red lines in Fig. 6), capturing the evolution from terrain-following (more stratiform) and less intense rainfall fields in the colder simulation cases to well-organized and intense bands for simulations at intermediate temperatures.For warmer simulations the rainfall pattern changes again to more intense but less organized (more cellular) convective features (Fig. 6d), and the UM parameters also show variation (Fig. 5c and d).These differences in convective dynamics have important effects in  3, and see also Kirhsbaum andDurran, 2005b andFuhrer andSchär, 2007).
The same analysis was performed on instantaneous horizontal cross sections of cloud water mixing ratio, q c , at two different levels (the height of q c maximum intensity, z q c max and z = 1100 m) after 3, 4 and 5 h of simulation, corresponding to the mature stage when the convective structures are fully grown.Figure 7 illustrates the robust scaling of moments of q c fields at z q c max for three different simulations.The same behavior is found for cloud fields in all simulations in Table 1, at all the analyzed times and heights.The UM parameters computed from q c also display complex nonlinear variation with simulation properties (red line at z q c max www.nonlin-processes-geophys.net/20/605/2013/ Nonlin.Processes Geophys., 20, 605-620, 2013  ac , red and blue lines are computed from q c at z q c max and z = 1100 m, respectively.Light blue markers are computed from ice mixing ratio fields at the level of their maximum intensity."x", "o", " " markers represent simulation times 05:00 h, 04:00 h and 03:00 h, respectively.and blue at z = 1100 m in Fig. 5).In particular the transition between colder and warmer cases is again clear.Both α and C 1 show similar values at the different times, due to the quasi-stationary state of the convective structures that is ac (red isolines), q c fields at z q c max at 04:00 h (black isolines) and q i at the level of its maximum intensity at 04:00 h (blue isolines) from simulations attained in these simulations.Comparison of the scaling parameters computed from q c fields at different vertical levels shows vertical variation of both UM parameters, although the shape of the curve is largely preserved.Notice that the z q c max level varies between simulations and with time, and it can sometimes capture different structures and cause significant variations of the parameters.This explains some values that seem to fall outside the curves: for example the value of C 1 Fig. 7. Relationship between M q and λ for different values of q.M q is computed from q c at z q c max at 04:00 h for simulation (a) Sst5km, (b) Sst10km and (c) Sst15km.The scales considered on the analysis are marked by "o" while "x" marks the ones not considered.The plots are on a log-log scale.
for Sst15km at 05:00 h computed from q c at z q c max is quite different from the one at 04:00 h, which is easily understood by examining Fig. 8: in Sst10km the z q c max level crosses the core of the cloud band, the same happening for Sst15km at 04:00 h, while at Sst15km at 05:00 h it is actually crossing a region of instabilities that has formed on the top region of the band, resulting in a different horizontal pattern and associated scaling parameters.The large variation in C 1 for the largest advective timescale τ adv at different times has a similar explanation.
The cumulative rainfall field analyzed here is viewed as representative of the entire simulated convective rain event.
It is important to stress that a time integrated quantity will depend on the evolution of the system during the integra- tion interval and the specific quantitative values of scaling parameters should change during this period due to the transient nature of the interactions among various physical processes.Presently, at this stage of the research, the analysis is focused on the spatial properties of the accumulated precipitation fields for a specific time of integration in a manner consistent with standard rainfall observations that consist of rainfall accumulations over a measurement timescale.Downscaling is critical to obtain robust spatial distributions of rainfall for use in hydrological studies and applications.

Scaling in cloud simulated fields -3-D analysis
The previous section showed a vertical variation of 2-D horizontal UM parameters.One possibility to avoid the complex problem of knowing the vertical variation of horizontal scaling exponents is to consider each horizontal level (different height) as a realization of the field and average over these several realizations, obtaining a single set of horizontal UM parameters representative of the horizontal scaling for the entire volume.The moments averaged over several vertical levels show improved linear fits as compared to the 2-D cases for all simulations (Fig. 9), and in particularly for higher order moments, which should be explained by the larger data sample size considered.Like in the 2-D case, some persistent deviations are identifiable but they have smaller amplitudes, and the particular scales for each simulation between 2-and 3-D cases are different, providing support to the notion of www.nonlin-processes-geophys.net/20/605/2013/ Nonlin.Processes Geophys., 20, 605-620, 2013 Fig. 9. Relationship between M q and λ for different values of q.M q is computed from q c at 04:00 h for simulations (a) Sst5km, (b) Sst10km, (c) Sst15km.The scales considered on the analysis are marked by "o" while "x" marks the ones not considered.The plots are on log-log scale.
absence of characteristic scales.The presence of scaling behavior for the entire volume implies that concurrent spatial structures (e.g., band core region, instabilities at the top of the bands, clouds that form downstream of the large scale orography, among others, Fig. 8) fall in the same average horizontal scaling regime.The scaling behavior in the vertical direction is not investigated here due to data limitations associated with the limited number of vertical levels in the model (see Sect. 2) together with the fact that we are looking at instantaneous realizations of fields rather than ensemble average statistics.But in other cases it could be investigated using a different set of UM parameters to describe anisotropy Green, blue and red lines are computed from q c .The light blue markers lines are computed from q i ."x", "o" and " " markers represent simulation times 05:00 h, 04:00 h and 03:00 h, respectively.between horizontal and vertical directions (see e.g., Lovejoy et al., 1987;Lazarev et al., 1994;Radkevich et al., 2007).
As in the 2-D case, the UM parameters vary with mean atmospheric and terrain properties (Fig. 10).Figure 10a shows C 1 decreasing with l sst while Fig. 10b shows there is little variation of α with l sst .Figure 10c and d show the abrupt transitions in both scaling parameters with changes in T s , further providing support for this result.Recall that in these idealized simulations the temperature profiles and cloud layer stability are determined by T s .However, in more realistic profiles this will not generally be true.Linear stability analysis suggests that a good alternative is to characterize the atmospheric environment using a (spatial) average N 2 m value.This average is computed at each time instant and only from points that exceed a fixed (arbitrary) q c threshold (e.g., 0.01 g kg −1 ) as it is defined only in saturated regions.Figure 11a and b show that the abrupt transition corresponds to the transition between the stable (N 2 m > 0) and unstable (N 2 m < 0) cases.For larger negative values, N 2 m < −2.5 × 10 −6 s −2 there is a break in variability of the UM parameters for small changes in N 2 m , particularly clear in C 1 .This behavior poses an intriguing m for all simulations.Black markers in plots (b) and (h) represent simulation Sst10km (i.e., the base case)."x", "o" and " " markers represent simulation times 05:00 h, 04:00 h and 03:00 h, respectively.question as to the character of the underlying physics and of this flow transition for highly unstable flows, not unlike the transition to high Reynolds numbers in the case of fully developed turbulence turbulence.The fields with the strongest instability values are not necessarily the warmer ones as defined by the initial conditions, as the average N 2 m value can be significantly altered by the release of latent heating or due to the presence of downstream stable clouds and other features.
To capture the effect of space-time co-localization of moist instability release and moist processes, N 2 m can be multiplied by a measure of the cloud fraction, cf, defined here by the ratio of number of grid points that exceed the chosen q c threshold to the total number of points.This transformation gives a measure of the overall realized moist instability at a given point in time, and provides a physically based framework to arrange the simulations in the order of colder (more stable, higher cf × N m 2 ) to warmer profiles (more unstable, lower cf × N m 2 ), consequently yielding more robust relationships (Fig. 11c and d).Indeed, this result holds for different values of the threshold, and the transition between stable and unstable cases becomes clearer in these plots.The stable region shows nearly constant UM parameters, with lower values of C 1 associated with a more space filling mean component in stratiform clouds as compared to more localized convective structures; the region of low realized instability corresponds to simulations that produce quasi-stationary rainbands, in agreement with the numerical experiments of Kirshbaum and Durran (2005a) who found that marginal potential instability and moderate wind speeds were the most favorable atmospheric conditions for convective bands to develop.This region displays complex variation of the scaling parameters that might be associated with the presence of different regimes due to competition between buoyancy and mechanical generation/dissipation of turbulence.Further investigation -with a larger population of simulations in this specific region of instability values -is necessary for further confirmation of this finding.In particular, the role of wind shear (usually understood as a mechanical turbulence generator) in band generation is not clear and requires more detailed analysis.The bands form in the absence of any basicstate vertical shear (as it is seen in our simulation and also in Kirshbaum and Durran, 2005a;and Kirhsbaum et al., 2007b), although wind shear will dynamically develop in these simulations eventually due, for example, to adjustment to terrain or formation of localized convective circulations.Kirshbaum and Durran (2005a) also found that in some cases, wind shear can actually suppress rainbands.Figure 11c and d also suggest that for higher instability conditions, when free convection dominates, and corresponding to the more cellular and disorganized patterns, the parameter variation becomes less complex.
The quantity cf × N m 2 is transient, depending on the dynamical evolution of the system.In particular, besides varying with upstream temperature profile, it also varies with τ adv and l sst .Thus we investigate its potential utility as a single predictor for the UM parameters for all the simulated q c fields analyzed here.Figure 11g and h show the UM parameters for all analyzed simulations.For α, all the points seem to fall into a single curve, suggesting cf × N m 2 might be used potentially as a single predictor.However C 1 shows variability with l sst (green and pink markers) and τ adv (red www.nonlin-processes-geophys.net/20/605/2013/ Nonlin.Processes Geophys., 20, 605-620, 2013 markers), particularly in the region of low instability (stable cases).Also the pink curves in Fig. 11c, d, and g and 11h show that, for the same initial temperature profiles but l sst = 5 km instead of 10 km, different values of C 1 are obtained.Thus cf × N m 2 alone cannot be used as a single predictor for both UM parameters, and variation with l sst and τ adv has to be considered separately.
As mentioned earlier, Perica and Foufoula-Georgiou (1996) found a linear dependency of statistical scaling parameters on another atmospheric stability measure, the Convective Available Potential Energy (CAPE), for a small number of rainfall fields of selected convective storms over relatively flat topography.This result was tested in our simulated fields which have a wider range of CAPE values, and a more complex nonlinear relation was obtained (Fig. 11e and  f).Compared to N 2 m , the use of CAPE has the disadvantage that the stable range of the profiles is ambiguous (CAPE = 0 for all stable cases), which presents a problem unless UM parameters were shown to be constant for all stable cases.Nevertheless, note that the unstable/stable transition of the parameters is still clear as CAPE goes to zero, as well as the complex nonlinear behavior for the small magnitude instability values -the banded cases -and again less complex behavior for higher magnitude instability values -cellular disorganized cases.This is consistent with Kirshbaum and Durran (2005b) who proposed that the distinction between more banded and cellular cases is related to enhanced susceptibility to free convection, with larger CAPE implying more cellular structures.They suggest that Convective INhibtion (CIN) could be an important parameter as well, but all our simulations have very small values of CIN and are not suitable for such analysis.

Relations between different scaling fields
In both 2-and 3-D analysis, UM parameters computed from P 5h ac and q c show clear relations for certain ranges of simulation properties, either in the shape of the variation curve in Fig. 5a, c, d and Fig. 10a, c, d, e, or for the exact values of the parameter, as can be seen for example in Fig. 5a for C 1 computed from q c at z q c max for l sst ≥ 8 km and C 1 computed from z = 1100 m for smaller l sst (also in Fig. 5c, see the C 1 computed from q c at z q c max for T s ≥ 285 K, and in Fig. 11c for T s ≥ 287.5 K).For the 2-D analysis these relationships depend on the chosen vertical level, sometimes best relating to the level of maximum intensity, sometimes to some other level, which makes its use difficult in practical applications, further supporting the 3-D analysis.
Figures 5c, d and 10c, d show that for colder simulations (T s ≤ 277.5 K) both UM parameters computed from P 5h ac are very similar to the ones computed from ice phase (snow + ice + graupel) mixing ratio fields, q i , which is particularly remarkable in the 3-D analysis.The role of cold microphysics on precipitation processes is also clear in simulation Sst10km_T272.5_WR(which is exactly the same as in Sst10km_T272.5,except that the Kessler warm rain scheme is used in WRF and ice microphysics is turned off).The simulation with warm rain processes only produces a cloud pattern similar to the one including cold microphysics, but it yields zero P 5h ac (Table 3).This result is strengthened by the fact that below 277.5 K the values of precipitation intensity and q i increase with decreasing temperature (Table 3).Finally, it is also apparent in Fig. 6a that for the colder case there is a clear relation between the P 5h ac (red line) and q i (blue line) patterns, while the q c (black line) seems to be more closely related to the underlying topography (grey scale).If the temperature is slightly raised, both q i and P 5h ac go to zero (Fig. 6b and Table 3).If the temperature is raised further, well-organized cloud bands form with a clearly associated P 5h ac pattern and no ice phase (Fig. 6c).Besides varying N 2 m , changing T s also changes the amounts of water in each phase and the vertical structure of the clouds, particularly the cloud depth (Fig. 12), which the linear theory also predicts as being an important parameter in embedded convection.For the warmer cases considered here, the cloud depth growths enough so that ice processes become important again (blue lines in Fig. 6d and Table 3), but the precipitation pattern now seems to be more related to q c (Fig. 6d).In these warmer cases, the parameter C 1 from P 5h ac seems to relate better to q c while α seems to relate better to q i (Fig. 10c, d).Overall, these results suggest that scaling of the simulated total precipitation fields is a complex function of the scaling of several other fields including water mixing ratio in the different phases, and probably also small-scale terrain spectra and τ adv .That is, the scaling of the integral precipitation fields reflects the underlying governing physical processes, which change with time and space as weather systems evolve and propagate over complex terrain.

Spectral analysis
The isotropic power spectra, described in Sect.4.1, can also be used to investigate the spatial scaling of horizontal fields.Here we average the spectra of instantaneous q c fields over several vertical levels, similar to what was done in Sect.4.3 for moment scaling analysis.Figure 13 shows three examples of log-log power spectra (black lines) computed at simulation hour 04:00: at the bottom for a colder stable case (Sst10km_T275), in the middle for an intermediate banded case (Sst10km), and at the top for a warmer more disorganized convective case (Sst10km_T292.5).Only the scales well resolved by the model are considered (and represented), assumed here to be larger than 5 x (limited by model resolution).It is important to keep in mind that Fourier analysis has the implicit requirement (assumption) of an infinite domain, and clearly here, as in general in all applications of Fourier analysis to digital data, we are limited by the domain size of the grid, to which the analysis is quite sensitive.It was found that the scale range we considered above L 0 /4 is too large based on the lack of symmetry of the domain, and the upper boundary of the range was reduced to a value L 0 /8.Outside this range, we do not have enough information on the power spectrum of the field.The black dashed lines in Fig. 13 are least-square fits to the considered range of scales (wavenumbers) of the power spectrum.There is a clear mean linear trend (scaling) in the log-log power spectra, with R 2 values for the least-square fits always above 0.90.Peaks in the power spectra, particularly evident in the intermediate banded case spectrum (middle spectra in Fig. 13), are consistent with harmonics resulting from the banded pattern.This behavior is expected when analyzing single realizations of cloud fields (e.g., Harris et al., 2001), and more so for events over complex terrain when orographic forcing is dominant (e.g., see scaling analysis in Barros et al., 2004).The scaling analysis was repeated here for all the simulations with l sst = 10 km and varying T s , all of them showing a mean linear trend in the log-log power spectra.Figure 14 shows that the variation of fitted values of β with cf × N 2 m displays similar transitions to the ones found for UM parameters, between stable/unstable cases and banded/disorganized convection.The location of the stable/unstable transition now seems to have been shifted to values slightly lower than N 2 m = 0.For numerical modeling parameterizations or downscaling, the exact location of this transition should be determined using a large data base of realistic simulations.Yet, the existence of such transition alone, which is supported by both spectra and moment scaling analysis and agrees with expected physical behavior such as for example the Richardson number criteria in boundary layer stability analysis, is an important finding in itself.The colder stratiform cases show some spread in β values that was not found in the UM parameters, including an unphysical β > 3 value.It should be taken into consideration these cases have shallower clouds (Fig. 12), which imply fewer Fig. 13.Log-log plots of isotropic power spectra computed from horizontal q c field at 04:00 h averaged over several vertical levels, from simulations: Sst10km_T292.5 (top), Sst10km (middle) and Sst10km_T272.5 (bottom).Only the range of scales between 5 x and about L 0 /8 is represented.The black dashed lines are least square fits to the spectra.levels of data for computation, along with low intensities of cloud fields, making them susceptible to computational uncertainty in the spectral exponent, adding to the limitations in the spectral analysis.The highly idealized topography with a superposed 2-D sinusoidal field used in the simulations here was designed to replicate and build upon previous works that used linear theory to investigate the development of orographic convection, and is not scaling.The structure of stratiform (terrain following) and orographically triggered convective fields is very influenced by the properties of the topographic elevation (e.g., Barros and Lettenmaier, 1994;Kirshbaum and Durran, 2004;Barros et al., 2004), and therefore a linkage between their scaling properties should exist.Follow up work, with realistic terrain should give further insight into this question.
Theoretically (isotropic) wind speed and passive scale advection in turbulent processes has β = 5/3 for 3-D turbulence, and 3 for 2-D turbulence.However cloud fields in convective turbulence should not be treated as passive scalars due to strong nonlinear effects such as latent heating release, entrainment, and active thermodynamics and microphysical processes, and therefore these theoretical results cannot be expected to hold.Additionally these theories assume isotropy between horizontal and vertical directions which is not the case in the atmosphere.

Summary and discussion
Linear stability analysis reveals that embedded convective structures are the result of unstable growth of small-scale disturbances, which should be very important in determining the predictability of orographic precipitation.They also suggest the governing role of parameters such as atmospheric stability, advective timescale, small-scale terrain spectra and cloud depth.The results from idealized numerical simulations of orographic convective precipitation qualitatively agree with these predictions.Particularly, they show that the representation of small-scale terrain (and associated lee waves) plays a very important role for these events, controlling the band spacing and significantly enhancing the convective intensity and rainfall amounts.But these simple models are unable to do any quantitative predictions or even predict the correct pattern (except in very particular cases) due to the presence of nonlinear effects, and alternative nonlinear models are required.
As pointed out by Schertzer and Lovejoy (1992), the contributions of small-scale activity to much larger scales are non-negligible, and indeed the corresponding variability has the most extreme behavior which explains in part the reason why multifractals have been successful on the study of extreme events.In fact several previous works reported statistical multiscaling behavior to be a very general property of geophysical fields (including rain and clouds), providing a convenient alternative framework due to its theoretical simplicity, its ability to handle nonlinear dynamics over very wide range of scales, its potential for downscaling applications and the fact that the convective events studied here are a particular kind of turbulent flow and hence a statisti-cal approach to the problem is appropriate.This multiscaling behavior holds for statistical moments and isotropic power spectra of the simulated rain and cloud fields in the considered range of scales (∼ 1 km to ∼ 20 km).The universal multifractal model reproduces quite well the empirical estimates of the horizontal scaling exponent function, K(q), and can be used to quantify the multiscaling behavior of the fields.Analysis of single horizontal sections of cloud fields revealed variation of the scaling parameters with vertical position.We argue that if we use an average over several horizontal sections at different vertical levels, we obtain estimates for the horizontal parameters of the entire volume, while the vertical structure should be represented by a different set of (vertical) scaling parameters.This transient 3-D anisotropic scaling framework seems to be a promising starting point to investigate relations between the statistical and physical properties and perform downscaling for mesoscale events.Results from both 2-and 3-D approaches revealed significant variations of K(q) and respective UM parameters with small-scale terrain spectra, atmospheric stability and advective timescale.Therefore, the representation of sub-grid scale variability of cloud and rainfall fields using statistical multifractal diagnostics would intrinsically depend not only on the properties of the field itself, but also on particular geographic location and atmospheric conditions, bringing the necessity of developing relationships to predict the scaling parameters from coarse grid atmospheric data and terrain spectra.However, the results also suggest that the development of such relationships is far from being trivial because they appear to involve complex nonlinear behavior.
A particularly robust behavior found here is the abrupt transition in the UM parameters between stable and unstable cases, which has a clear physical correspondence as a transition from a more stratiform to a more organized (banded) convective regime.This behavior is also found in the spectral exponent β.Another instability measure, CAPE, also seems to be a good predictor, but it is unable to distinguish between the different stable cases and CIN might have to be used in conjunction.N 2 m also has the advantage of allowing a linkage with the analytical models where it appears explicitly, which is not so clear for CAPE.It was also found that the scaling of the surface accumulated precipitation field (usually an important field to be parameterized and downscaled) is a complex function of the horizontal scaling of several other fields such as water mixing ratios in the different phases and probably l sst and τ adv .
The scaling of the power spectrum can be taken advantage to perform stochastic spatial downscaling (disaggregation) of instantaneous rain and cloud field (see e.g., Bindlish and Barros, 2002;or Rebora et al., 2006).The present work suggests that scaling parameters could be computed from known large scale measures (such as stability and terrain spectra).Although the focus of spectral analysis is on second order statistics, and hence the scaling behavior at other orders is neglected by this type of downscaling approach, the time Nonlin.Processes Geophys., 20, 605-620, 2013 www.nonlin-processes-geophys.net/20/605/2013/ dependency of the scaling parameters is introduced by the temporal evolution of the atmospheric conditions, and the usually small time step used in numerical simulations can be advantageous to achieve a high temporal representation.Further improvements could include scaling constraints using higher-order moments.Finally, it is important to point out some limitations of our study.In order to investigate the dependence of the multifractal properties on mean properties, a limited range of highly idealized simulations are used where we can control the mean atmosphere with a small number of simulation parameters.Because our work suggests a transient dependence of the scaling on mean properties, the analysis is performed on instantaneous fields, rather than ensemble averages, which increases the volume of data necessary for analysis.Horizontal isotropy is assumed, and it is important to stress that the issue of existence of horizontal anisotropy must be tackled in the future.The solutions show robust horizontal scaling even though the only external forcing is a very simple terrain with poor scaling and the initial conditions are horizontally homogeneous.In realistic cases, the atmosphere and terrain forcing will themselves be scaling, which might influence the scaling of the solutions.Models such as WRF involve many subgrid scale parameterizations which vary among model configurations.Here we use parameterizations only for subgrid turbulence closure and microphysical processes.Whether these parameterizations introduce unphysical scaling behavior in the models at resolved scales, and whether the scaling depends on the particular chosen parameterizations are legitimate questions.Furthermore, in the real atmosphere there are many other effects not considered here such as radiation, earth rotation, surface and boundary layer effects which might introduce further complications and influence the scaling.In this study, however we spell out unambiguously what processes are resolved or not, and there is general a body of literature that shows that the resolution adopted can resolve explicitly dominant convective structures (see for example Fuhrer and Schär, 2007).In fact, as pointed out in Sect.2, our simulations are based on previous work where the authors showed that they reproduce the essential features of the observed convective structures in the Oregon Coastal Range.Harris et al. (2001) compared the scaling observed and model produced precipitation fields at 3-km resolution, finding an encouragingly similar behavior at scales above 5 times the resolution (i.e., 15 km).Although they simulated realistic cases with different model and parameterizations than here, the model is based on the same governing equations.Our assumption here is that our simulations, based on nonlinear, compressible and non-hydrostatic Navier-Stokes equations coupled with the thermodynamics equations in the WRF model under the specified forcing conditions represent the key physical processes of interest and that the multifractal parameters can be used to describe the behavior of the simulated system in the explicitly well-resolved scales.It is important to stress that the fields resulting from numerical simulations should reproduce the scaling in geophysical observations if they are to be realistic.Furthermore, we apply multifractal analysis not as means to find universal parameters that will exhibit consistency, though that could be a result, but rather as a mathematical model can capture, at any given scale, variability that reflects a large dynamical range of the phenomena determined by the long-range dependencies.Thus, what is consistent is the model used to describe variability; we then make use of the multifractal model parameters as metrics that can help interpret model dynamics and thermodynamics in a systematic way.The presence of an underlying multifractal behavior is robustly supported by the evidence presented: a clear mean linear trend in the spectral analysis in conjunction with the robust scaling in the moment analysis, which is less sensitive to the data limitations and reveals robust scaling behavior up to orders of about 2, as it is typically found in atmospheric fields and reported in literature consistent with the scope of the present manuscript.
Although much insight was gained into scaling dependencies, particularly the presence of transient dependence of scaling parameters on mean conditions found here, an investigation on observational atmospheric data and realistic numerical simulations is warranted, along with further support to whether numerical weather prediction models are able to reproduce, at least to some extent, the multifractal properties of simulated fields and what is the effect of the particular parameterizations chosen on the scaling.This is currently ongoing work.

Fig. 1 .
Fig. 1.Horizontal maps of topographic elevation (in meters) for simulations: (a) CTL and (b) Sst10km.The inner black square represents the domain of interest where the scaling analysis computations are performed.
(ii) the large-scale orography is a finite length ridge and open boundary conditions are used in this study for all lateral boundaries instead of a quasi-1-D ridge with periodic boundary conditions in the y boundaries.Periodic boundary conditions used in previous studies would reintroduce small-scale www.nonlin-processes-geophys.net/20
Fig.3.Relationship between M q and λ for different values of q.M q is computed from P 5h ac from simulation (a) Sst5km, (b) Sst10km, (c) Sst15km.The scales considered on the analysis are marked by "o" while "x" marks the ones not considered.The plots are on a log-log scale.

Fig. 4 .
Fig. 4. The scaling exponent function, K(q) for (a) varying l sst , (b) varying T s with l sst = 10 km and (c) varying U with l sst = 10 km.The empirical estimated K(q) is showed with markers and the UM model fit is represented by full line.

Fig. 5 .
Fig. 5. Relationship between UM parameters (C 1 and α) and: (a, b) l sst ; (c, d) T s ; and (e, f) τ adv .Black lines represent values computed from P 5hac , red and blue lines are computed from q c at z q c max and z = 1100 m, respectively.Light blue markers are computed from ice mixing ratio fields at the level of their maximum intensity."x", "o", " " markers represent simulation times 05:00 h, 04:00 h and 03:00 h, respectively.

Fig. 8 .
Fig. 8. Cross sections (xz) of q c at the y position of its maximum intensity for simulations (a) Sst10km at 04:00 h, (b) Sst15km at 04:00 h and (c) Sst15km at 05:00 h.The solid black line represents topography.The upper and lower black dashed lines represent the levels z q c max and z = 1100 m, respectively.

Fig. 10 .
Fig. 10.C 1 and α computed from 3-D averages of q c horizontal fluctuations against: (a, b) l sst , (c, d) T s and (e, f) τ adv .Black lines represent values computed from P 5hac .Green, blue and red lines are computed from q c .The light blue markers lines are computed from q i ."x", "o" and " " markers represent simulation times 05:00 h, 04:00 h and 03:00 h, respectively.

Fig. 11 .
Fig. 11.C 1 and α computed from 3-D averages of q c horizontal fluctuations against: (a, b) N 2 m ; (c, d) cf × N 2 m ; and (e, f) CAPE.Blue lines and markers are computed from simulations with varying T s and l sst = 10 km, pink from simulations with varying T s and l sst = 5 km, red for varying τ adv and green for varying l sst .Plots (g, h) show C 1 and α against cf × N 2m for all simulations.Black markers in plots (b) and (h) represent simulation Sst10km (i.e., the base case)."x", "o" and " " markers represent simulation times 05:00 h, 04:00 h and 03:00 h, respectively.

Fig. 12 .
Fig. 12. Vertical profiles of q c at the position of their maximum intensity at 04:00 h, for simulations with l sst = 10 km and varying T s .

Table 1 .
Summary of WRF simulations population.

Table 2 .
Local maximum and domain total of P 5h ac field in simulations with different small-scale terrain wavelengths.

Table 3 .
Local maximum of P 5h ac , q c and q i fields and domain total P 5h ac for simulations with different temperature profiles.