D Koutroumpa

R Lallement

J C Raymond

V Kharchenko

THE SOLAR WIND CHARGE-TRANSFER X-RAY EMISSION IN THE 1/4 keV ENERGY RANGE: INFERENCES ON LOCAL BUBBLE HOT GAS AT LOW Z

Keywords: interplanetary medium -ISM, bubbles -ISM, general -supernova remnants -X-rays, diffuse background -X-rays, general -X-rays, ISM

We present calculations of the heliospheric solar wind charge-exchange (SWCX) emission spectra and the resulting contributions of this diffuse background in the ROSAT 1/4 keV bands. We compare our results with the soft X-ray background (SXRB) emission detected in front of 378 identified shadowing regions during the ROSAT All-Sky Survey. This foreground component is principally attributed to the hot gas of the so-called Local Bubble (LB), an irregularly shaped cavity of ∼50-150 pc around the Sun, which is supposed to contain ∼10 6 K plasma. Our results suggest that the SWCX emission from the heliosphere is bright enough to account for most of the foreground emission toward the majority of low galactic latitude directions, where the LB is the least extended. On the other hand, in a large part of directions with galactic latitude above 30 • , the heliospheric SWCX intensity is significantly smaller than the measured one. However, the SWCX R2/R1 band ratio differs slightly from the data in the galactic center direction, and more significantly in the galactic anticentre direction where the observed ratio is the smallest. Assuming that both SWCX and hot gas emission are present and their relative contributions vary with direction, we tested a series of thermal plasma spectra for temperatures ranging from 10 5 to 10 6.5 K and searched for a combination of SWCX spectra and thermal emission matching the observed intensities and band ratios, while simultaneously being compatible with O vi emission measurements. In the frame of collisional equilibrium models and for solar abundances, the range we derive for hot gas temperature and emission measure cannot reproduce the Wisconsin C/B band ratio. This implies that accounting for SWCX contamination does not remove these known disagreements between data and classical hot gas models. We emphasize the need for additional atomic data, describing consistently EUV and X-ray photon spectra of the charge-exchange emission of heavier solar wind ions.

INTRODUCTION

The diffuse soft X-ray background (SXRB), first observed in the 1970s (Bowyer et al. 1968;Williamson et al. 1974;Sanders et al. 1977) has since been shown to be the sum of local and distant sources. Above 2 keV, it is dominated by the extragalactic background, itself a combination of unresolved point sources and warm-hot interstellar medium (WHIM) diffuse emission (Hasinger et al. 1993). At lower energies, it is dominated by the galactic halo (Burrows & Mendenhall 1991;Snowden et al. 1994), and finally below 0.3 keV it is mainly due to the unabsorbed emission from hot gas filling, the so-called Local Bubble (LB; McCammon et al. 1983;Bloch et al. 1986;Snowden et al. 1990aSnowden et al. , 1990b)), a cavity devoid of dense gas extended at high latitudes and connected to the halo (Frisch & York 1983;Welsh et al. 1998;Lallement et al. 2003). The main tools used to disentangle local and distant emission are the "shadowing" experiments, i.e., spatial variations of intensity and spectral characteristics around and toward dense, soft X-ray absorbing clouds (e.g., Herbstmeier et al. 1995). Snowden et al. (1998Snowden et al. ( , 2000) ) used more than 370 ROSAT shadows to produce almost full-sky mapping of the "unabsorbed" component of the emission, i.e., the LB contribution.

This was the generally accepted scenario until the discovery of X-ray emission in comets (Lisse et al. 1996) and the identification of the emission mechanism as charge-exchange (CX) reactions between the highly charged heavy solar wind (SW) ions and the cometary neutrals (Cravens 1997). Cox (1998) suggested that the CX reactions should also occur between heavy SW ions and interstellar neutrals (H and He) in interplanetary space and that the resulting X-ray emission (solar wind charge exchange; SWCX) should have an impact on the SXRB interpretation. Cravens (2000) estimated that the quiescent level of SWCX emission could be of the same order as the SXRB component attributed to the LB. This interplanetary, heliospheric emission is time-dependent because of the intrinsic variable nature of the SW. Short timescale variations tend to be washed out by integration along the line of sight (LOS) and the size of the emitting region (Cravens et al. 2001), but longer term variations, including those related to the solar cycle, can cause more persistent changes in the heliospheric emission level. In addition, there is a contribution from the Earth's magnetosphere, due to charge transfer with exospheric neutrals. Such emission, studied in detail by Robertson et al. (2006), reacts instantaneously to SW variations and magnetosphere shape variations, leading to a high variability and the occurrence of high intensity peaks following solar events. The spectral characteristics of some spectacular enhancements have been recently recorded by the XMM-Newton and Suzaku satellites (Snowden et al. 2004;Henley & Shelton 2008). For such events both magnetospheric and heliospheric contributions may be present.

Very likely most of the sharp increases of terrestrial origin have been removed from the ROSAT map along with the cleaning procedure of the Long Term Enhancements (LTEs Snowden et al. 1994), as well as some heliospheric increases, especially toward the downwind (DW) side of the interstellar flow where the gravitational cone of focused helium is the most reactive region. Indeed, most points on the sky in the ROSAT map were observed several times over the course of at least two days, allowing identification and removal of periods of enhanced emission. A debate is still maintained, though, about the actual level of the quasi-stationary heliospheric contribution to the ROSAT maps of unabsorbed emission. This contribution is extremely difficult to detect from time variations. On the other hand, SCWX and hot gas thermal emission have different spectral properties, i.e., the observed spectral information should help to disentangle the two processes. This is the subject of the present study. For a recent review of all types of SWCX phenomena see Bhardwaj et al. (2007).

The first estimates of the stationary heliospheric contribution (Cravens 2000;[START_REF] Lallement | [END_REF]) were based on simplifying assumptions about the spectral characteristics of the SWCX emission. Since then the existence of the SWCX phenomenon has motivated theoretical work on exact photon yield values for the charge-transfer collisions (Kharchenko & Dalgarno 2000;Pepino et al. 2004), as well as a number of laboratory experiments devoted to the CX emission mechanism. For a recent review see Wargelin et al. (2008).

The spatial distribution of magnetospheric and heliospheric SWCX emission was modeled by Robertson & Cravens (2003), revealing significant variations in brightness as a function of Earth location, LOS direction, and activity phase. Koutroumpa et al. (2006) computed similar maps for a few specific energy bands after including Pepino et al.'s (2004) detailed CX emission spectra for C, N, O, and Ne ions. [START_REF] Lallement | [END_REF], taking into account the specific viewing geometry of ROSAT showed that the heliospheric background in the 1/4 keV band was nearly isotropic and could have been unnoticed in the All-Sky Survey maps, while accounting for a large portion of the signal, and possibly the major part at low galactic latitudes.

Using both the stationary and time-dependent models, Koutroumpa et al. (2007Koutroumpa et al. ( , 2009) ) modeled four high-latitude shadowing observations and showed that in the 3/4 keV band, where the oxygen lines (O vii triplet at 0.57 keV and O viii line at 0.65 keV) are dominant, the SWCX emission from the heliosphere can account for all the unabsorbed, local component of the SXRB, with no need of an LB emission. In parallel, the SW contribution to the background and its variability have been shown to be responsible for some discrepant measurements (Smith et al. 2007) and for supposedly low-energy counterparts of distant objects (Bregman & Lloyd-Davies 2006).

A 10 6 K plasma, however, has very little emission in the 3/4 keV band and mainly emits in the 1/4 keV band. The Koutroumpa et al. (2007) results, while not requiring any LB emission, therefore do not preclude the existence of 10 6 K gas. Exact calculations of the SWCX spectra and intensities below 0.3 keV are mandatory if one wants to disentangle LB hot gas diffuse emission from the SWCX background. In this paper, we examine the SWCX contribution to the 1/4 keV spectral region, compare this contribution to observations and also to contributions from hot gas at different temperatures.

Independently of the SWCX contribution, a number of results have somewhat contradicted the interpretation of the unabsorbed SXRB as the LB 10 6 K gas emission.

1. Data from the NASA Extreme Ultraviolet Explorer (EUVE) satellite and from the dedicated CHIPS mission did not detect the EUV emission expected from surrounding 10 6 K gas (Jelinsky et al. 1995;Hurwitz et al. 2005). It has been suggested that a very low metal abundance may be responsible for this nondetection, but the required depletion level corresponds to the physical state of very dense clouds, which is unlikely for 10 6 K, tenuous gas. 2. The pressure of this hot gas derived from the X-ray background is far above the pressure within the local interstellar cloud and other clouds embedded in the LB [START_REF] Lallement | The Local Bubble and Beyond[END_REF]Jenkins 2002). 3. Low-latitude absorption measurements of highly charged ions such as Si iv, C iv, and O vi formed in conductive interfaces between the hot (10 6 K) gas and embedded cold:warm clouds do not seem to correspond to expectations from the models (Slavin & Frisch 2002;Indebetouw & Shull 2004). Column densities of Si iv and C iv are too small and line widths too narrow (Welsh & Lallement 2005), and O vi is detected only at the periphery of the local cavity, while one would also expect interfaces between the hot gas and the local clouds (Welsh & Lallement 2008). 4. Fundamental discrepancies arise also when comparing the Wisconsin sounding rocket survey data in the B and C bands, and the ROSAT All-Sky Survey data in the R1 and R2 bands. The four bands are pictured in Figure 1 (upper panel). In the low energy (0.1-0.2 keV) B band (Bloch et al. 1986;Snowden et al. 1994), the intensity seems to be higher than what is predicted by thermal emission models. This has been particularly well demonstrated by [START_REF] Bellm | [END_REF] who have made a global study over the 0.1-0.3 keV interval. According to this work, a best fit to all energy bands is provided by very low metallicity gas at 10 5.85 K, but inspection of their results (see their Figure 4) reveals significant discrepancies between measured and observed ratios for this best-fit solution.

Especially, the B/R12 band ratio favors a low temperature 10 5.8 K (the B band intensity is high and favors a shift of the spectrum toward low energies), while the R2/R1 ratio favors temperatures above 10 6 K (R2 is relatively high, favoring a shift toward high energies).

Whether or not the existence of the SWCX background can help to explain part or all these contradictions is a question that has now to be addressed. This work is a first step in this direction. In Section 2, we describe the SWCX emission and spectral model we have developed and how we make use in our analysis of the Raymond and Smith (R-S) hot plasma model. In Section 3, we compute the expected SWCX emission and the contribution in ROSAT R1 and R2 bands for each of the 378 shadow regions observed by Snowden et al. (2000). We compare the SWCX R1 + R2 intensity with the unabsorbed component derived by the Snowden et al. (2000) shadow analysis and discuss the distribution of the discrepancies between data and the SWCX model. In Section 4, we compute the SWCX model R2/R1 and B/C band ratios, as well as the corresponding ratios for hot gas (R-S model) in collisional equilibrium within a large temperature range. We compare the modeled ratios with the observed band ratios during the two (Wisconsin and ROSAT) surveys. In Section 5, we search for a combination of SWCX and hot gas emission compatible with the observed intensities and band ratios and we compare those solutions with observational constraints from O vi and EUV background measurements. In Section 6, we discuss the results and draw some conclusions. ). The emitting ions are marked above the most prominent lines.

MODEL DESCRIPTION

SWCX Model

The basic model calculating the SWCX emission in the inner heliosphere was thoroughly presented in Koutroumpa et al. (2006); parameters appropriate for the ROSAT All-Sky Survey are discussed in Section 3. We calculate self-consistently the neutral H and He density distributions in the inner heliosphere (up to ∼100 AU), in response to solar gravity, radiation pressure, and anisotropic ionization processes for the two neutral species. Ionization of H atoms is mainly due to their CX collisions with SW protons and He atoms are mostly ionized by solar EUV photons and electron impact. We also consider the impact of CX on the SW ion distributions. This interaction is described in the following reaction:

X Q+ + [H, He] → X * (Q-1)+ + [H + , He + ].
(1)

The collision rate per volume unit R X Q+ (cm -3 s -1 ) of X Q+ ions with the neutral heliospheric atoms is given by the equation:

R X Q+ (r) = N X Q+ (r) υ r (σ (H,X Q+ ) n H (r) + σ (He,X Q+ ) n He (r)) = R (X Q+ ,H) (r) + R (X Q+ ,He) (r), (2)
where σ (H,X Q+ ) and σ (He,X Q+ ) are the hydrogen and helium CX cross sections, n H (r) and n He (r) are the hydrogen and helium density distributions, respectively, ῡr = VSW -ῡn ≈ VSW the relative velocity between SW ions and interstellar (IS) neutrals in the inner heliosphere, and N X Q+ (r) is the self-consistent solution to the differential equation:

dN X Q+ dx = -N X Q+ (σ (H,X Q+ ) n H (x) + σ (He,X Q+ ) n He (x)) + N X (Q+1)+ (σ (H,X (Q+1)+ ) n H (x) + σ (He,X (Q+1)+ ) n He (x)) (3)
expressing the evolution of the density distribution of ion X Q+ along SW streamlines due to production (from CX reactions of ion X (Q+1)+ ) and loss terms.

Cross section uncertainties are mainly due to instrumental systematic errors and most important to the collision energy dependence of cross sections. Detailed uncertainties for individual ions are not given in literature, but average uncertainties of ∼30% at most are reported (Wargelin et al. 2008).

Then, we establish emissivity grids in units of (photons cm -3 s -1 ):

ε i (r) = R (X Q+ ,H) (r) Y (E i ,H) + R (X Q+ ,He) (r) Y (E i ,He) , ( 4
)
where Y (E i ,M) is the photon emission yield (in number of photons) computed for a spectral line of photon energy E i following CX with the corresponding neutral species M (H or He individually). For any LOS and observation date, the directional intensity of this spectral line is given by

I E i (LU ) = 1 4π ∼100 AU 0 ε i (s) ds, (5)
which defines the average intensity, in line units (LU = photons cm -2 s -1 sr -1 ), of the spectral line for the particular date and LOS, as well as the solar cycle phase (minimum or maximum) corresponding at this date. The intensity is somewhat underestimated because of the SW ion propagation in the heliosheath up to the heliopause, and in the heliotail up to ∼3000 AU, where all ions are used up. The outer heliospheric region is neglected in our model, but estimates yield a maximum additional ∼20% contribution in the DW direction, with possible effects on the SWCX spectral hardness (see Section 4).

Our original atomic database (Kharchenko 2005) included C 5,6+ , N 5,6,7+ , O 6,7,8+ , Ne 8,9+ , and Mg 10,11+ ions. Exact calculations of the cascading photon spectra were performed individually for these ions when they CX with hydrogen and helium, respectively. Detailed CX collision cross sections taking into account both the neutral target species and the SW velocity regime were include in the calculations (P. Stancil 2003, private communication). These calculations have already been used to reproduce observed SWCX spectra from comets with CHIPS (Sasseen et al. 2006).

The database was recently updated to include Fe 7...13+ , Si 5...10+ , S 6...11+ , and Mg 4...9+ ions that emit intense lines in the 0.1-0.3 keV range. Individual emission spectra induced in the CX collisions of these ions have very complicated structures because of a large number of intermediate multiplets related to different excited states of many-electron ions. The photon yields Y(E,M) for heavier ions were calculated using the simplified hydrogenic model (Kharchenko & Dalgarno 2000), which assumes a hydrogenic nature of electronic excited states. In this model, the effective charge of hydrogenic ion is computed from an accurate value of the ion ionization potential and branching ratios of radiative cascading transitions are chosen to be the same as in all H-like ions. Moreover, photon yields were calculated using a single neutral species, which means that no distinction between H and He was made. The hydrogenic approximation of the CX emission spectra is a quantum mechanical model in which an actual ion spectra may be replaced with the hydrogenic spectra. In this model, the total energy of emitted photons is defined by an initial state population and should be an accurate quantity matching real spectra. Positions of emission lines do not correspond exactly to real emission spectra, but this defect is not very important at the low resolution of the observed spectra. Total cross sections of CX collisions for the hydrogenic approximations have been calculated using the overbarrier model (Kharchenko & Dalgarno 2001). An example of calculated spectra is presented in Figure 1(lower panel), with the emitting ion identifying the most intense lines.

Hot Gas Thermal Emission

We use an R-S hot plasma model (Raymond & Smith 1977) assuming typical metal abundances [He, C, N, O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni] = [10. 93, 8.52, 7.96, 8.82, 7.96, 7.52, 7.60, 7.20, 6.90, 6.30, 7.60, 6.30] [START_REF] Allen | Astrophysical Quantities[END_REF]. We use this model rather than the APEC model (Smith et al. 2001) that superseded it because we are most interested in the 1/4 keV range. APEC includes only transitions for which accurate atomic rates are available, while the code of R-S estimates the emission in the large number of weak lines that are known to be present (e.g., from moderately ionized species of Mg, Si, S, and Fe) but which lack accurate excitation rates and wavelengths. Given the low spectral resolution of the observations considered here and our interest in the total emitted power, the R-S model serves very well.

This model gives us X-ray emissivities f 1 (T) and f 2 (T) convolved by and summed in the ROSAT R1 and R2 bands, respectively, as a function of temperature such that the total hot gas X-ray flux in these bands is defined as

I 12,LB = I 1,LB + I 2,LB = EM(i, T ) • (f 1 (T ) + f 2 (T )), (6)
where EM(i, T) is the emission measure for temperature T and look direction i. Units of functions f 1 and f 2 are RU EM -1 , where RU = 10 -6 counts s -1 arcmin -2 is the usual ROSAT detector unit and EM is the typical emission measure unit cm -6 pc.

Equivalently, the hot gas emission in bands B and C is defined as

I B C,LB = I B,LB + I C,LB = EM(i, T ) • (f B (T ) + f C (T )), (7)
where f B (T) and f C (T) are the equivalent emissivity functions in the B and C bands derived by the R-S plasma code, in units of counts s -1 EM -1 .

R1 + R2 INTENSITIES

We have calculated SWCX spectra in ROSAT observation geometry for the shadow field LOS listed in Table 1 of Snowden et al. (2000), which were observed during the ROSAT All-Sky Survey. The shadows analyzed by Snowden et al. (2000) were located at high galactic latitudes and in general above 15 • from the galactic plane. The ROSAT observation geometry is defined with the view direction perpendicular to the Sunsatellite direction. Thus, it takes a six-month period to build a full-sky map of the soft X-ray intensity. The ROSAT All-Sky Survey was performed between 1990 July and 1991 February, which corresponds to maximum solar activity conditions that were taken into account in the SWCX simulations.

Maximum solar activity conditions imply the following input parameters in the SWCX model. We consider a radiation pressure to gravity ratio μ = 1.5 for neutral hydrogen and slightly anisotropic ionization rates varying between 8.4 × 10 -7 s -1 at the solar equator and 6.7 × 10 -7 s -1 at the poles (Quémerais et al. 2006). For neutral helium, the average lifetime (inverse ionization rate) at 1 AU is 0.62 × 10 7 s at solar maximum, in agreement with McMullin et al. (2004). In solar maximum, the SW is considered to be a complex mix of slow and fast wind states, that is, in general approximated with an average slow wind flux. Slow SW flows at ∼400 km s -1 and has a proton density of ∼6.5 cm -3 at the Earth position. The oxygen content with respect to protons is [O/H] = 1/1780. The most important heavy ion charge state abundances with respect to oxygen [X q+ / O] are C 5,6+ :[0.21, 0.318], O 6,7,8+ :[0.73, 0.2, 0.07], Si 8,9,10+ :[0.057, 0.049, 0.021], and Fe 8,9,10,11+ :[0.034, 0.041, 0.031, 0.023] (adopted from Schwadron & Cravens 2000).

We have convolved the individual spectra with the ROSAT R1 and R2 band responses in order to calculate the total SWCX flux in these bands, as well as the total R12 (R1 + R2) flux. We plot the resulting R12 SWCX flux and the unabsorbed I 12,obs component from the Snowden et al. (2000) analysis as a function of absolute galactic latitude in Figure 2. I 12,obs corresponds to the unabsorbed portion of the SXRB that was originally attributed to the LB ∼ 10 6 K hot gas. X-ray intensities are presented in ROSAT units (RU = 10 -6 counts s -1 arcmin -2 ).

The SWCX R12 flux (black dots) varies between 212 and 460 RU with an average value of 332 RU and is fairly uniformly distributed across all latitudes. The lower and upper limits calculated in the SWCX simulations for average maximum conditions are represented with the plain black horizontal lines.

On the other hand, the unabsorbed I 12,obs component (gray circles) derived in the Snowden et al. (2000) analysis has a clear correlation with the absolute galactic latitude. Higher I 12,obs values are measured toward higher latitudes, where the local cavity is enlarged and communicates with the galactic halo through the chimneys.

In the figure, it is clear that the SWCX intensity is of the same order as the I 12,obs intensity measured in low galactic latitudes (up to around 20 • -25 • ). We can conclude, then, that the SWCX 1/4 keV flux could account for most of the observed ROSAT emission in the galactic plane. This conclusion implicitly assumes that the highly peaked exospheric SWCX contribution has been cleaned from the ROSAT data, but not the heliospheric contribution.

Figure 3 is a map in galactic coordinates of the ratio between the unabsorbed ROSAT emission and the computed SWCX contribution. The map clearly reveals the emission from the so-called chimneys that connect the local cavity to the northern 2000), gray circles) and calculated SWCX (black dots) fluxes in the ROSAT 1/4 keV (R1 + R2) band as a function of absolute galactic latitude for 378 shadowing observations. Flux is given in ROSAT Units (10 -6 counts s -1 arcmin -2 ). and the southern halo. Our intensity results do not preclude that outside these chimneys the totality of the signal is SWCX emission.

BAND RATIOS

For each SWCX spectrum calculated in the look directions presented in Figure 2, we have calculated the R2/R1 (ROSAT) and C/B3 (Wisconsin) ratios. We find an average R2/R1 (hereafter R CX ) ratio of 1.39 and an average C/B ratio of 6.67. Although both ratios show very little variation across the sky, there is a hardness trend of the SWCX spectra with harder spectra toward the DW direction (UW to DW variations: R CX = [1.36-1.41], C/B = [6.25-7.14]). However, we need to alert the reader that these are somewhat uncertain SWCX spectra and therefore somewhat uncertain band ratios. Indeed, as we mentioned in Section 2, exact calculations for Fe, Si, Mg, and Al

are not yet available, and no distinction was done between the neutral targets (H or He), while laboratory experiments show that the energy levels populated after the electron capture and the subsequent radiative cascades may differ significantly for different targets. Since most of the SWCX DW emission is due to the interaction with neutral helium, while on the upwind (UW) side hydrogen is the main contributor, more precise calculations could have an effect on the hardness. Also, although preliminary calculations show an almost negligible effect, the heliospheric model cutoff (especially in the DW directions) may be responsible for the "loss" of relatively more emission from lower charge states (at relatively lower energies) than emission from higher charge states (at relatively higher energies). Thus, the calculated SWCX spectra may actually be softer than what is predicted here. It is evident that a more detailed calculation taking into account all metals and the neutral target nature, as well as detailed cascading collisions (secondary ion production) in the outer heliosphere is needed in the future. On the other hand, the interval we find for the ratio can be used a reliable value for the average SWCX ratio.

These SWCX ratios have to be compared with the corresponding ratios for thermal emission. The latter were obtained as a function of temperature by convolving the R-S spectra with the ROSAT band responses R1 and R2 and Wisconsin B and C responses. The results are shown in Figure 4. Above log T = 6.1 the thermal R2/R1 ratio reaches its maximal value of 1.2. It remains, however, slightly lower than the SWCX ratio of 1.36-1.41. At those temperatures, the thermal C/B ratio increases to its maximal value of 4, a value almost half the SWCX ratio of 6.24-7.14, i.e., a significant difference. Those curves allow estimates of the ratios for combinations of thermal plus SWCX background emissions.

COMBINATION OF THE HELIOSPHERIC SWCX AND LB HOT PLASMA EMISSION

For a comparison with the data we consider two regions: one centered on the direction of the incoming IS flow at ecliptic coordinates (λ, β) ∼ (252.3, 8.5) • for the IS H flow, according to Lallement et al. (2005) outgoing flow direction (DW direction). In galactic coordinates, the UW direction corresponds to (l, b) = (5.1, 19.6) • , close to the galactic center direction (antigalactic direction for DW, respectively). These two regions are also very close to the minimum and maximum values of the hardness ratio derived by Snowden et al. (1990b), which define the so-called color gradient axis of the SXRB. For these two regions, we can derive average values of observed unabsorbed 1/4 keV emission using the LB contours in the Snowden et al. (1998) analysis. For the UW direction, the ROSAT unabsorbed I 12,obs emission we estimate ∼325 RU, while for the DW direction the observed unabsorbed level is found to be 450 RU.

The equivalent B+C intensities in the Wisconsin survey are estimated on average ∼90 and ∼125 counts s -1 for the UW and DW directions, respectively (Snowden et al. 1990b). However, those intensities include both the foreground (assumed LB) and more distant components (galactic halo and extragalactic), since the Wisconsin survey did not have enough spatial resolution to study the shadowing fields. Moreover, the Wisconsin sounding rocket measurements looked in the roughly anti-Sunward direction, which should also affect the comparison with the ROSAT All-Sky Survey in terms of the SWCX component spatial distribution. For instance, for the DW look directions, the Wisconsin sounding rockets were observing directly through the He cone and should have had a higher "contamination" of SWCX emission than the ROSAT detectors that must have been located in crosswind positions on the Earth's orbit in order to observe in the DW directions.

The measured R2/R1 and C/B ratios for these regions are shown superimposed to the models in Figure 4 and listed in Table 1. For the UW area, the measured R2/R1 ratio is close to the SWCX value, although slightly smaller, while for the DW area it is significantly smaller. For both areas, the measured C/B ratio is lower than the SWCX ratio. Figure 4 shows that for both areas, a combination of SWCX and thermal emissions may, in principle, account for the ROSAT measurements, the thermal emission lowering the R2/R1 ratio to achieve the observed value. Similarly, independently of the R2/R1 ratio, Figure 4 also shows that a combination of both backgrounds may account for the Wisconsin data, the thermal emission decreasing the C/B ratio to achieve the observed value.

It remains to find a combination satisfying both ratios simultaneously. Our attempt to find a solution is the following one. For each assumed temperature of the hot gas, we use the R1 and R2 data (and thus the observed ratio) and the SWCX spectral shape to derive the respective contributions of SWCX and thermal emission, i.e., we derive which quantity of SWCX induced R12 intensity and which emission measure EM for the hot gas led to the measured intensities and the measured R2/R1 ratio. We a posteriori calculate the B and C intensities, C/B ratio, and the O vi emission of the hot gas and compare with the data.

The SWCX model predicts a total R12 intensity:

I 12,CX = I 1,CX + I 2,CX = I 1,CX • (1 + R CX ) ( 8 ) = I 2,CX • 1 + R CX R CX , ( 9
)
where R CX is the R2/R1 ratio predicted by the SWCX model. The total unabsorbed flux I 12,obs (i) measured in the R12 band toward look direction i is the sum of LB hot gas I 12,LB and SWCX I 12,CX fluxes: I 12,obs (i) = I 12,CX + I 12,LB , so that SWCX intensity can be written as I 12,CX = I 12,obs (i) -I 12,LB = I 12,obs (i) -EM • (f 1 (T ) + f 2 (T )).

(10) The observed R2/R1 ratio (hereafter R obs ) toward look direction i is defined by the equation

R obs (i) = I 2,CX + I 2,LB I 1,CX + I 1,LB . ( 11
)
Resolving Equation (11) by using Equations ( 6)-( 10), we find the hot gas emission measure EM(i, T) as a function of temperature, total R12 measured intensity I 12,obs (i) and measured R obs (i) ratio toward look direction i:

EM(i, T ) = R CX -R obs (i) (R obs (i) + 1) • (R CX • f 1 (T ) -f 2 (T ) )
•I 12,obs (i).

(12) We calculate the emission measure for the two UW and DW directions defined above using the following numerical values:

(1) the R CX ratio is constant and equal to 1.39, (2) observed values of the unabsorbed portion of the 1/4 keV emission in the R12 band are I 12,obs (UW, DW) = (325, 450) RU as derived from the LB contours in the Snowden et al. (1998) analysis for the UW and DW (respectively galactic and antigalactic) directions, and (3) the corresponding observed R2/R1 ratio is R obs = 1.25 and 1.04 for the UW and DW directions, respectively.

In Figure 5, we show the resulting EM and the portion of the total emission due to the SWCX mechanism (called residual emission) as derived from Equations ( 12) and (10), respectively, as a function of log T and for the two look lines. We show superimposed the R12 measured intensities in those directions.

For the calculated EM and corresponding temperatures we have also added to Figure 5 the intensity of the O vi doublet at λ = 1034 Å (1032 and 1038 Å). In order to calculate this O vi doublet emission, we have used Equation ( 5) of Shull & Slavin (1994) and assumed that interstellar O abundance is 8.5 × 10 -4 . We also assume that the O vi ion proportion depends on temperature according to the Chianti database formulas for collisional equilibrium (Landi et al. 2006).

In order to delimit the possible temperature solutions for the LB hot gas, we place the following constraints: (1) We assume that the SWCX model is accurate enough to ensure that the heliospheric emission in the R12 band cannot be lower than ∼212 RU (lower limit in Figure 2). This gives us (from the right panel of Figure 5) an upper limit of log T = 6 in temperature.

(2) We use the observed upper limit of O vi doublet intensity, reported at ∼800 LU (Shelton 2003), which gives us two temperature limits at log T = 5.13 and 5.64 (extreme limits in the DW direction). The interval between those two temperatures is forbidden because the corresponding O vi column densities (and thus the O vi intensity) are too high to match observations. Temperatures below 10 5.13 K would predict extremely strong C vi and N vi absorption toward nearby stars, which have not been observed (e.g., Lehner et al. 2003) so we do not consider it a realistic solution. The limits of valid temperature intervals are marked by the vertical bold lines and the cross-hatched regions in Figure 5 show the excluded temperature ranges. The two most plausible hot gas temperature limits (log T = 5.64, 6.0) along with the corresponding UW and DW emission measures and residual SWCX emission are summarized in Table 1.

We also calculate for the two boundary solutions the corresponding SWCX intensities and the thermal emission intensities in the B and C bands by convolving our simulated SWCX and the hot gas spectra with the band responses. For the two temperatures, the total hot gas and SWCX intensity in the B + C band (I BC in Table 1) is found to be about 83 counts s -1 and about 115 counts s -1 for the UW and DW directions, respectively. This similarity arises from the similarity between the wavelength intervals covered by the B and C bands and the R1 and R2 bands (see Figure 1). The C/B ratio derived from this analysis is 5.0 and ∼3.5 for the UW and DW directions accordingly, for temperatures above 10 5.64 K.

The (B + C) total intensity is consistent with the lower values reported in the Wisconsin survey (Snowden et al. 1990b), which correspond to the lower galactic latitudes. Moreover, Snowden et al. (1990b) did not proceed with a shadowing analysis of the Wisconsin data, so the reported values include both local and more distant absorbed components and are expected to be higher than the hot gas and SWCX combination we present here.

However, the C/B ratio computed in the analysis (5 to ∼3.6 from UW to DW, depending on temperature) is inconsistent with the observed value, especially in the DW direction (observed ∼2.2), suggesting that we should need more hot gas emitting in the B band. This inconsistency cannot be attributed to the absorbed portion of emission included in the Wisconsin data analysis because the absorbed component is a high-T gas giving a harder spectrum since absorption is more effective in lower energies.

This inconsistency of the DW C/B ratio is important, since it seems difficult to explain in the context of our study. As a matter of fact, as can be seen in Table 1, the SWCX contribution in the C band is large whatever the temperature, and reaching a C/B ratio of 2.2 requires a very small SWCX emission, in our sense far from realistic. Again, as we discussed in the introduction and it was shown in the [START_REF] Bellm | [END_REF] study, the B intensity is higher than expected from the models. This seems to remain true (and even worse) when taking into account the SWCX contribution.

DISCUSSION AND CONCLUSIONS

We have modeled the intensity and spectral characteristics of the heliospheric SWCX emission at the time of the ROSAT survey and compared with the unabsorbed, local emission derived by Snowden et al. (2000) in the 1/4 keV band. The results show that the SWCX emission can account for most of the total intensity recorded in the R1 + R2 bands for most of low-latitude LOS. A map of the heliospheric SWCX portion of the total signal clearly reveals the high-latitude chimneys to the halo as the only regions unambiguously dominated by hot gas emission. Such a result can be interpreted as meaning that little or no hot gas exists within the galactic disk.

The spectral characteristics, however, reveal more complexity and preclude such a simple scenario. The SWCX band ratios disagree with the observations, especially toward the galactic anticenter and at low energies (C/B). We have thus searched for a combination of SWCX and thermal emission from hot gas in equilibrium and solar abundances able to reproduce the data. Our study shows that a combination of SWCX and thermal emission can reproduce the data in the galactic center hemisphere at low latitudes. For this solution, the SWCX emission strongly dominates. The temperature of the hot gas is constrained within the interval 10 5.64 -10 6 . The upper limit is constrained by the lower limit on the SWCX intensity. This upper limit can be considered as firmly determined, thanks to recent observational studies above 0.3 keV that have confirmed the validity of our model (Koutroumpa et al. 2007(Koutroumpa et al. , 2009)). The temperatures lower than 10 5.64 are excluded by O vi emission observations (Shelton 2003) and interstellar ion absorption lines toward nearby stars.

On the other hand, it is difficult to fit with such a combination the Wisconsin data. In the UW (galactic center direction), a combination of hot gas and SWCX emission gives (B + C) intensities as well as C/B ratios roughly compatible with the observed values for several different temperature ranges. The main difficulty is the impossibility to account for the very low C/B ratio measured toward the galactic anticenter direction with the present input models used in our study. We note that the high B intensity is also clearly a problem for any hot gas solution, including the very high depletion hypothesis, as shown by the study of [START_REF] Bellm | [END_REF]. The SWCX contribution, which hardens the spectra, reinforces this difficulty.

However, further investigation is required on the model's uncertainties in order to quantify their influence on the spectral hardness of SWCX emission. Further analysis of the hydrogenic ion approximation is needed, since hydrogenic ions tend to emit photons at higher energies following CX (since high-n to ground transitions are generally allowed) than the multielectron ions considered here (because selection rules and more complicated atomic structure lead to more cascades before the final transition to ground). Therefore, it is likely that the model spectra have significantly less flux at low energies than they should. Given the steeply decreasing effective areas at lower energies (see Figure 1), this would help explain some of the discrepancies seen in the C/B ratio.

Also, the fact that no distinction was made between the neutral targets (H or He) in the 1/4 keV calculations, would also have an effect on the spectral hardness, since there is an effect on the electron capture level (roughly proportional to I n -1/2 , where I n is the neutral target ionization potential), the capture level with He being somewhat lower than with H. As noted in Section 5, the Wisconsin survey is more sensitive to this effect than ROSAT because of its observation geometry (i.e., the Wisconsin DW observations looked more directly through the He cone, where SWCX emission has a softer spectrum), again helping to explain some of the C/B ratio anomaly. Finally, at the low energies of the B band, more detailed calculations of the heliospheric and magnetospheric signals must be performed, especially the low-energy secondary SWCX emissions in the heliosheath and heliotail, i.e., subsequent recombinations of partially neutralized SW high ions.

NoFigure 1 .

 1 Figure1. Top panel: Wisconsin B (plain) and C (dashed), and ROSAT R1 (dash-dotted) and R2 (dotted) band effective areas. Bottom panel: example of calculated SWCX spectra in line units (photons cm -2 s -1 sr -1 ). The emitting ions are marked above the most prominent lines.

NoFigure 2 .

 2 Figure 2. Measured unabsorbed (I 12,obs , fromSnowden et al. (2000), gray circles) and calculated SWCX (black dots) fluxes in the ROSAT 1/4 keV (R1 + R2) band as a function of absolute galactic latitude for 378 shadowing observations. Flux is given in ROSAT Units (10 -6 counts s -1 arcmin -2 ).

Figure 3 .

 3 Figure 3. Unabsorbed I 12,obs R12 flux data over simulated I 12,CX SWCX flux ratio map in galactic coordinates with an equal-area Aitoff projection. For information, we note the UW and DW directions in the sky. The Local Cavity (LC) chimneys appear clearly where the I 12,obs /I 12,CX ratio is the highest. (A color version of this figure is available in the online journal.)

Notes.

  In the lower part of the table, we include the observational input considered.a I B = EM • f B (T), I C = EM • f C (T). b I 12,CX = I 12,obs -I 12,LB . c I BC = I B (LB + CX) + I C (LB + CX), C/B = I B (LB + CX) / I C (LB + CX).

Figure 4 .

 4 Figure 4. Wisconsin (upper panel) and ROSAT (lower panel) band ratios for data (UW bold, DW plain line), SWCX model (dashed line) and the R-S hot gas model (dotted line) as a function of temperature.

NoFigure 5 .

 5 Figure 5. Solutions for the EM (dashed line) and residual SWCX emission (dot-dashed line) as a function of temperature of a hot plasma for observed R12 intensities (I 12,obs , plain horizontal lines) for an upwind (left panel) and downwind (right panel) look directions. Upper O vi doublet intensity limit (dotted line) is calculated for the EM values. Discontinuities in the O vi intensity curve are due to O vi density interpolation in the Chianti database. The cross-hatched regions show the hot gas temperature ranges excluded by O vi intensity (Shelton 2003) and SWCX emission constraints.

Table 1

 1 Temperature (log T) Limits and EM Solutions for Upwind (UW) and Downwind (DW) Look Directions When Combining an SWCX and R-S Hot Plasma Code

	Look		Local Bubble							SWCX		(LB + SWCX) c
	Direction	log T	E M	I 12,LB	I B	a	I C	a	I 12,CX	b	I B	I C	I BC	C/B
			(10 -4 cm -6 pc)	(RU)	(counts s -1 )	(RU)	(counts s -1 )	(counts s -1 )
	UW	5.64	5.5	25	3.8	2.6	300		10.	66.5	82.9	5.
	DW	21.	96	14.5	10.	354	11.8		78.8	115.1	3.37
	UW	6.00	3.7	62	4.8	11.1	263		9.2	58.	82.7	5.
	DW		14.1	238	18.3	42.5	212		7.7	46.5	114.1	3.57
			Observational Input
		I 12,obs (RU)	R2/R1 (R obs )	I BC (counts s -1 )	C/B
	UW	325	1.25	∼ 90	4.
	DW	450	1.04	∼ 125	2.17

The original papers on the Wisconsin survey referred to the B/C ratio, but given the extensive use of the R2/R1 ratio in our analysis and the rough correspondence between R1 and B, and R2 and C, we refer to the C/B ratio.

The authors wish to thank Dan McCammon and Steve Snowden for comments and useful discussions on the ROSAT and Wisconsin surveys and band responses. The work presented has benefited from very useful discussions at the ISSI workshop "From the Outer Heliosphere to the Local Bubble" and at the "Local Bubble and Beyond II" meeting. We are extremely grateful to the anonymous referee for his attentive report and valuable comments which resulted in significantly improving the paper.


Stefano Federico

Chiara Bertolin
email: chiara.bertolin@ntnu.no

Nazzareno Diodato

Pasquale Borrelli

Panos Panagos

Gianni Bellocchi

Communicating Hydrological Hazard-Prone Areas in Italy With Geospatial Probability Maps

Keywords: erosive density, hydrological hazard, prone area, geospatial probability map, Italy

The recurrence of storm aggressiveness and the associated erosivity density are detrimental hydrological features for soil conservation and planning. The present work illustrates for the first time downscaled spatial pattern probabilities of erosive density to identify damaging hydrological hazard-prone areas in Italy. The hydrological hazard was estimated from the erosivity density exceeded the threshold of 3 MJ ha -1 h -1 at 219 rain gauges in Italy during the three most erosive months of the year, from August to October. To this end, a lognormal kriging (LNPK) provided a soft description of the erosivity density in terms of exceedance probabilities at a spatial resolution of 10 km, which is a way to mitigate the uncertainties associated with the spatial classification of damaging hydrological hazards. Hazard-prone areas cover 65% of the Italian territory in the month of August, followed by September and October with 50 and 30% of the territory, respectively. The geospatial probability maps elaborated with this method achieved an improved spatial forecast, which may contribute to better land-use planning and civil protection both in Italy and potentially in Europe.

INTRODUCTION

Changes in the frequency of severe weather events determine the occurrence of damaging hydrological events [START_REF] Diodato | A millennium-long reconstruction of damaging hydrological events across Italy[END_REF][START_REF] Petrucci | MEFF: the database of Mediterranean flood fatalities (1980 to 2015)[END_REF], which in turn affect many natural and human systems [START_REF] Corella | Three millennia of heavy rainfalls in Western Mediterranean: frequency, seasonality and atmospheric drivers[END_REF][START_REF] Sofia | Flood dynamics in urbanised landscapes: 100 years of climate and humans interaction[END_REF]. Modeling rainfall intensity can help to quantifying the effect of these extreme phenomena that involve complex landscape and ecosystem scenarios within multiple event feedback forms (e.g., [START_REF] Thomas | Landscape sensitivity in time and spacean introduction[END_REF][START_REF] Mulligan | Modelling and model building[END_REF][START_REF] Harris | Biological responses to the press and pulse of climate trends and extreme events[END_REF].

Water is certainly a resource for the Earth's ecosystems, but intense rainfall and overland flows can also be land disturbing forces. Extreme hydro-meteorological events such as rainstorms produce high-impact land processes, which increases with their intensity [START_REF] Toy | Soil Erosion: Prediction, Measurement, and Control[END_REF][START_REF] Wei | Effects of rainfall change on water erosion processes in terrestrial ecosystems: a review[END_REF]. However, diverse parts of the terrestrial ecosystems respond in different ways to unequal forcing agents, often non-linear and universal, with threshold-type characteristics [START_REF] Thomas | Landscape sensitivity in time and spacean introduction[END_REF][START_REF] Arnell | Uncertainty in the relationship between climate forcing and hydrological response in UK catchments[END_REF][START_REF] Knight | The impacts of climate change on terrestrial earth surface systems[END_REF][START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF]. Simultaneous mapping of these thresholds, and the processes for developing them, is a challenge in geographical information sciences [START_REF] Srivastava | Threshold concepts in geographical information systems: a step towards conceptual understanding[END_REF][START_REF] Li | Spatio-temporal variation of wet and dry spell characteristics of tropical precipitation in Singapore and its association with ENSO[END_REF]. The scale and the relative uncertainty integration of these thresholds, for which storms drive overland flows, including the soil erosion by water, are important because they affect the spatial pattern of damaging hydrological events, and ultimately improve our mechanistic understanding of ecological responses to climate extremes [START_REF] Kayler | Experiments to confront the environmental extremes of climate change[END_REF]. An example of rainfall thresholddriven process is the initiation of landslides [START_REF] Guzzetti | Rainfall thresholds for the initiation of landslides in central and southern Europe[END_REF]. Concepts such as thresholds of change and magnitude are all important since landscapes may be able to resist or absorb impulses of change as a form of sensitivity or, inversely, stability [START_REF] Huggett | Dissipative systems: implications for geomorphology[END_REF][START_REF] Meyer | Quantifying resilience to recurrent ecosystem disturbances using flow-kick dynamics[END_REF]. The application of this concept to the way how geomorphological systems respond to climatic variability and change is nowadays a relevant topic [START_REF] Allison | The sensitivity of landscapes[END_REF][START_REF] Phillips | Evolutionary geomorphology: thresholds and nonlinearity in landform response to environmental change[END_REF]. In geomorphology, the first use of sensitivity has been in the form of a ratio between resisting to and disturbing forces, inherently changing, but natural systems generally have a high capacity to absorb such changes [START_REF] Brunsden | A critical assessment of the sensitivity concept in geomorphology[END_REF]. Research indicates that distinct ways in which sensitivity analysis is applied have different implications, because whereas the ratio of resisting to disturbing forces is dimensionless, landforms within a landscape may react differently to environmental events, possibly owing to varied proximity to a particular threshold, which introduces an element of probability related to the likelihood of change [START_REF] Downs | The sensitivity of river channels in the landscape system[END_REF]. Thus, in this work the erosivity density (MJ ha -1 h -1 mm -1 ), that is, the ratio of rainfall erosivity to precipitation, is introduced to measure the erosivity per rainfall unit. The critical value of erosivity density ≥3 MJ ha -1 h -1 has been introduced in storm geomorphology by [START_REF] Dabney | Enhancing RUSLE to include runoff-driven phenomena[END_REF]. This hydrological threshold could help identify erosionand overland flow-prone areas [START_REF] Kinnell | Event soil loss, runoff and the universal soil loss equation family of models: a review[END_REF][START_REF] Panagos | Monthly rainfall erosivity: conversion factors for different time resolutions and regional assessments[END_REF]. The advantages of using erosivity density over directly calculated rainfall erosivity are (after [START_REF] Renard | Universal soil loss equation[END_REF]: (a) higher susceptibility to classifying hydrological hazards; (b) higher stability, achieved with a shorter period of records; and (c) easier interpolation over large territories owing to its altitudeindependency up to 3,000 m a.s.l. [START_REF] Foster | Comparison of the USLE, RUSLE1. 06c, and RUSLE2 for application to highly disturbed lands[END_REF].

Uncertainty in erosivity density thresholds poses challenges for the spatial analysis since the worst storms may fall in areas between recording weather stations (e.g., [START_REF] Willmott | Rising estimates of terrestrial and global precipitation[END_REF]. The summer season in Mediterranean lands is characterized by aggressive storms developing the risk for damaging hydrology and catastrophic events [START_REF] Pereira | Modelling extreme precipitation in hazardous mountain areas: contribution to landscape planning and environmental management[END_REF][START_REF] Diodato | How the aggressiveness of rainfalls in the Mediterranean lands is enhanced by climate change[END_REF][START_REF] Diakakis | Public perception of flood risk in flash flood prone areas of Eastern Mediterranean: the case of Attica Region in Greece[END_REF][START_REF] Salvati | Gender, age and circumstances analysis of flood and landslide fatalities in Italy[END_REF]. Local interactions between environmental context and climate variability are provided as comprehensive explanation for landscape changes and impact evaluation [START_REF] Bintliff | Time, process and catastrophism in the study of Mediterranean alluvial history: a review[END_REF]. The diverse geomorphology characterizing the Mediterranean region has important consequences for the atmospheric and sea circulations, with an uneven distribution of weather types [START_REF] Lionello | Cyclones in the Mediterranean region: climatology and effects on the environment[END_REF] and a broad spectrum of hydro-geomorphological hazard potentials across Italy [START_REF] Diodato | A millennium-long reconstruction of damaging hydrological events across Italy[END_REF].

Throughout the Mediterranean basin (including Italy), rainfall is typically concentrated between autumn and spring [START_REF] Trigo | Relations between variability in the Mediterranean region and mid-latitude variability[END_REF] but extreme hydrological and geomorphological events occur mainly in periods of imbalance of precipitation input, from May to October. Though only a few studies have explicitly investigated how precipitation affects the damaging hydrological events in this vulnerable period of the year [START_REF] Cevasco | Storminess and geo-hydrological events affecting small coastal basins in a terraced Mediterranean environment[END_REF], in recent years, a considerable progress has been made in developing rainfall erosivity databases from regional [START_REF] Borrelli | Rainfall erosivity in Italy: a national scale spatio-temporal assessment[END_REF] to continental scale [START_REF] Panagos | Rainfall erosivity in Europe[END_REF] and even to global [START_REF] Panagos | Global rainfall erosivity assessment based on high-temporal resolution rainfall records[END_REF] scale. However, the development of regional-scale spatial models of the storm hazard remains a challenge that scientists face for hydrological, agricultural, and environmental planning. The spatial uncertainty associated with the danger of extreme rainfall across a range of scales is an open issue, which is added to uncertainties associated with downscaling methods and lack of high temporal resolution data in many regions.

The present paper has the main objective of demonstrating to which extent the threshold associated with erosivity densities depends on hydrological damaging events prediction of. A coupled GIS and geostatistics model is the tool, applied here, to continuously map the spatial uncertainty of erosivity density thresholds and predict prone areas to damaging hydrological events in Italy. The achieved results are funded on climatic data and on a probabilistic approach. Other parameters, affecting the intensity and the frequency of damage in risk-prone areas as population density, infrastructures, plants density, and lowintensity rainy periods have not considered as out of the scopus of the present work.

MATERIALS AND METHODS

Study Area and Data Collection

The Mediterranean has the highest amounts of precipitation during the period October to March, while mean summer precipitation is less predictable owing to high rainfall variance in this season [START_REF] Xoplaki | Wet season Mediterranean precipitation variability: influence of large-scale dynamics and trends[END_REF]. Moreover, the period between end of summer and beginning of autumn results in considerable percentage of bare and tilled soils (compared to green cover) across the Mediterranean, including Italy. This agro-environmental condition accelerates the hazard of severe erosion and overland flow [START_REF] Pereira | Modelling extreme precipitation in hazardous mountain areas: contribution to landscape planning and environmental management[END_REF]. In this period, the sparser the plant cover, the more vulnerable the topsoil is to both dislodgement and removal by raindrop impact and surface runoff. As with spatial variability, the temporal variation in precipitation and vegetation response are nonlinearly related, highlighting the existence of threshold behaviors, which are heterogeneous in dynamic, changing landscapes [START_REF] Snyder | Multi-scale temporal variation in water availability: implications for vegetation dynamics in arid and semiarid ecosystems[END_REF].

The Italian territory is characterized by heterogeneous geographic systems, with sea surface waters around the Italian peninsula, where vast mountain chains, coexist with valleys and plains. The most southern part of Italy is an active storm area, and its orography is an important agent in modifying the track of cyclones and their synoptic evolution across Europe [START_REF] Homar | Sensitivities of an intense Mediterranean cyclone: analysis and validation[END_REF]. Such heterogeneity modifies mesoscale components of the circulation and generates precipitation patterns fluctuating greatly in space or time [START_REF] Luterbacher | A review of 2000 years of paleoclimatic evidence in the Mediterranean[END_REF]. Convective rainfall events, in particular, cause extreme, and erosive storms, which are geographically localized in Italy. Such large spatial and temporal variability can be difficult to capture accurately with rain gauge networks or remote sensing observations [START_REF] Curić | Comparison of modeled and observed accumulated convective precipitation in mountainous and flat land areas[END_REF]. Advective rainfall, in contrast, is characterized by weak rainfall intensity over large areas.

In this study, data on rainfall erosivity and precipitation were derived from the Rainfall Erosivity Database on the European Scale (REDES) [START_REF] Panagos | Rainfall erosivity in Europe[END_REF]. 219 uniformly distributed meteorological stations were selected in order to cover the different biogeographical and climate zones of Italy (Figures 1A,B) with an average density of about one station per 1,000 km 2 : The stations' elevation ranges from 0 to 2,541 m a.s.l. with a distribution reflecting the altitudinal gradients of the country. The frequency of acquisition of precipitation data is hourly over a period larger than 10 years (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). The gridbased map of rainfall erosivity density was created interpolating the data resources represented by 219 stations.

Erosivity Density Model

In an erosive climate, the magnitude, intensity, and frequency of a given precipitation has an impact in damaging hydrological processes as it changes both landform and the equilibrium of environmental systems. The assumption, in order to model the erosivity density (ED), is that the environmental system adapts to changes of the natural hydrological regime. In fact a fluctuation in such regime, especially when thresholds for an acceptable level of hydrological disturbance are exceeded, may have damaging consequences to the environment.

The ED algorithm, which predicts damaging hydrological hazards, is based on the ratio between storm erosivity and precipitation [START_REF] Dabney | Enhancing RUSLE to include runoff-driven phenomena[END_REF][START_REF] Renard | Universal soil loss equation[END_REF]:

ED = 1 n n i=1 m j k=1 I 30 • m r=1 0.29 • 1 -0.72e (-0.05i r ) • v r • k f • 1 n n i=1 P j (1)
where n is the number of recorded years; m j is the number of erosive events during a given month j; k is the index of a k th single event; v r is the rainfall volume (mm) during the r th period of a storm, which splits into m parts; I 30 is the maximum 30-min rainfall intensity (mm h -1 ); i r is the rainfall intensity during the time interval (mm h -1 ); P j is the rainfall amount (mm) during a given month j.

Therefore, the numerator in Equation ( 1) represents the mean monthly storm erosivity (MJ mm ha -1 h -1 ), while the denominator the relative average precipitation amount (mm) multiplied by a factor f.

The variables involved can be subjected to large time fluctuations and set threshold in the natural hydrogeomorphological regime, in such a case the coefficient f is a useful parameter to explain the degree of local ecosystem flexibility, whose estimate is required for the geomorphological risk assessment within given features [START_REF] Diodato | Geostatistical uncertainty modelling for the environmental hazard assessment during single erosive rainstorm events[END_REF]. In fact, natural land-based ecosystems can be quite flexible and capable of absorbing the stresses provoked by disturbance in a vast array of processes [START_REF] Mendoza | Integrating multicriteria analysis and GIS for land condition assessment: part I -evaluation and restoration of military training[END_REF], including the damages associated with weather extremes [START_REF] De Luís | Climatic trends, disturbances and short-term vegetation dynamics in a Mediterranean shrubland[END_REF], so that f >1. Contrarily, highly disturbed or degraded landscapes (e.g., intensive cropland, chaotic urbanization, post-fire vegetation dynamics, sparse vegetation), are commonly less flexible, so that f < 1.

According to [START_REF] Dabney | Enhancing RUSLE to include runoff-driven phenomena[END_REF], when ED>3 MJ ha -1 h -1 , storm erosivity accounts for a large percentage of the rainfall amount in an intense event, highlighting an increasing erosive hazard (above the critical value for damaging hydrological events).

Probability kriging and Fit of the Regional Model

Geostatistical techniques consider that distant data values are less related to each other than near spatial data values (after Tobler, 1970). The commonly used ordinary kriging (OK), in particular, computes an unsampled value (z), knowing its coordinates (x, y), and neighbor values. OK makes optimal, unbiased estimates of regionalized variables at unsampled points from values of the same variables at surrounding stations, using structural analysis and the initial set of measured data [START_REF] Journel | Mining Geostatistics[END_REF][START_REF] Journel | Nonparametric estimation of spatial distributions[END_REF]. Kriging-type approaches are probabilistic in nature because the uncertainty associated with a local estimate can be calculated. OK of a set of indicatortransformed values provides a resultant value between 0 and 1 for each point estimate. This is an estimate of the proportion of the values in the neighborhood, which are greater than the indicator or threshold value. Probability kriging (PK) was proposed by [START_REF] Sullivan | Nonparametric estimation of spatial distributions[END_REF] to overcome some issues raised with indicator kriging (IK) which is a non-parametric approach to estimate data uncertainty and is used to predict conditional probability (e.g., probability of the occurrence of an event above/below a threshold) for categorical data of unsampled location. Probability kriging is a non-linear method, which uses the information from all the available indicator variables by using the order relation of observed variables. It represents an attempt to obtain estimates less sensitive than the ones obtained with IK to the choice and number of cut-off values (thresholds). In addition, compared to IK, the estimates of PK fully reflect local variability [START_REF] Hohn | Geostatistics and Petroleum Geology[END_REF]. PK can be used to classify areas under damaging hydrological hazard prone areas. A straightforward approach is to classify as hazardous all locations where the probability of exceeding a critical threshold value, z k , is greater than a critical value of erosivity density (ED = 3 MJ ha -1 h -1 ). Kriging methods rely on the notion of autocorrelation. Correlation is usually thought of as the tendency for two or more variables to be related. For example, in climate sciences, weather stations close to each other have similar climatic values. When weather stations become distant, then climatic values are different. In the same way, in geostatistics, locations distant from measured climatic data may have high prediction uncertainty. The rate at which the correlation decays can be expressed as a function of distance and at some distance to weather stations, the pixel values have no correlation with the prediction location. Based on this assumption, a search neighborhood function (here below cdf) is generally specified (see result section) that limits the number and the configuration of the measured data points that will be used in the prediction. The mechanism controlling the radius of the neighborhood function shape is mainly driven from the range, which assumes different shapes in different months. In this way, PK uses more information than IK because the rank of the datum z(u α ) in the sample cumulative distribution function (cdf) is taken into account in addition to the indicator value exeeceding the threshold value z k [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF] as in the simplest IK approach.

Subsequently, the trade-off cost for this conditioning of the posterior cdf is the inference of the covariance function C X (h) and the K cross-covariances CC XI (h; z k ) of the model of coregionalization between indicator and uniform transforms at each threshold value (h is the distance from the station u α ). The coregionalization is the covariance matrices describing the correlation structure of the set of erosivity density indicator variable at a characteristic spatial scale. Cross-covariances give the covariances of variables at pairs of points.

The coindicator transform i(u;z k ), using the rank-order transform x(u) as a secondary variable, is the probability kriging estimator [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF]:

[I (u; z k )] * = n(u) α=1 λ α (u; z k ) • I (u α ; z k ) + n(u) α=1 υ α (u; z k ) • X (u α ) (2)
where the weights λ α (u; z k ) and υ α (u; z k ) are determined by solving the (2n(u) + 2) ordinary cokriging system of equations referred in [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF]. Since erosivity density data presents a non-Gaussian distribution, a logarithmic transformation was performed to reduce skewness before applying Equation (2). Then, logarithms were used for the probability kriging, named lognormal probability kriging (LNPK). The Geostatistical Wizard routine used to perform the spatial elaboration is in [START_REF] Johnston | Using ArcGis Geostatistical Analyst[END_REF] and runs under ESRI ArcMap (http://desktop.arcgis. com/fr/arcmap).

In addition, to fit the data to a regional model, a twostage iterative procedure developed by [START_REF] Johnston | Using ArcGis Geostatistical Analyst[END_REF] Finally cross-validation and error assessment procedures check for the internal consistency of the model and the spatial correspondence of damaging hydrological events have been carried out. The error involved during the interpolation of point data to landscape values through probability kriging is assessed through a quantitative estimation of the standard error and of the cross-validation.

RESULTS AND DISCUSSION

Structural Analysis and Modeling

In the Italian peninsula the summer erosivity density (ED) is, on average, 2.5 times larger than for the rest of the year [START_REF] Borrelli | Rainfall erosivity in Italy: a national scale spatio-temporal assessment[END_REF][START_REF] Panagos | Monthly rainfall erosivity: conversion factors for different time resolutions and regional assessments[END_REF]. The rainstorm are therefore considerably higher than the rest of Europe as the mean erosivity density in Italy is 2.02 MJ ha -1 h -1 , compared to 0.92 MJ ha -1 h -1 in Europe. In addition, although Northern Italy is particularly stormy, extreme rain events are also common in Central and Southern Italy [START_REF] Fiorillo | Reconstruction of a storm map and new approach in the definition of categories of the extreme rainfall, northeastern Sicily[END_REF]. This manifest the fact that, as we move from south to the driest hydro-climatic areas, the storm rainfall intensities remain high. This trend is reflected by the mean annual precipitation against the ratio of the maximum recorded daily to mean annual precipitation (Figure 1C). In fact, the highest ratios correspond to the driest locations. For instance, in areas where the mean annual rainfall is below 600 mm this ratio is typically around 10%, that is, an important fraction of the total annual rainfall which fall during very few storms. As a result, where rain is less abundant, the flood risk may become higher than in wetter situations where people are more familiar with frequent rainfall.

Based on the specified neighborhood function cdf (Figure 2), the months of August and September result in similar ranges i.e., a similar mechanism or driving factor controlling the radius of cdf shape, both reaching a search neighborhood with semi-axis of about 50 km (Figures 2A,B). The range of October is instead larger, reaching a search neighborhood with semi-axis of about 100 km (Figure 2C).

Trend analysis shows that only a weak gradient of the erosivity density data occurs along the northwest to southeast direction. The non-existence of a non-random (deterministic) component can be attributed to the persistent rainfall, which is received from the Atlantic. Based on this finding, the stationarity hypothesis cannot be rejected both locally and for the whole country.

Then, the implementation of the iterative procedure provides a satisfactory fit of covariance and cross-covariance models (Figure 3). This fit results from summing the geostatistical nugget and an exponential model with the previous estimated ranges of 50 km, for August and September, and 100 km, for October (solid curve in Figure 3, respectively). The obtained covariance and cross-covariance functions have been then used to train the geospatial tool on how control points act during the interpolation with LNPK. According to the isotropic neighborhood of data-points, five points per target-site have been used for estimation, corresponding to the search circle of the erosivity density-threshold (z k =3 MJ ha -1 h -1 ).

More specifically, this exponential model assumes that, at the monthly scale, the range variations of the erosivity density threshold are associated with the storm characteristics and local geographical features. This suggests that the increasing variability of hydrological events produce a more regular pattern of spatial variation in October and a spot-mosaic pattern in August and September. These months suggest that an enhanced hazard is associated with complex, more localized storm events, characterized by a strong convective component. In addition, the nugget ratio is larger for the August month (0.073), than to the months of September and October (around 0.05). This represents the random (unexplained) variance caused by the inherent variability in the data, and the spatial variation at distance much shorter than sample spacing which cannot be resolved in August at the scale where sampling is performed (< 20 km), or owing to some measurement errors.

Spatial Patterns of Erosivity Density Exceeding Threshold Values

As already reported in [START_REF] Ballabio | Mapping monthly rainfall erosivity in Europe[END_REF], the most erosive months in Italy are August (in North Italy), September (in Central Italy) and October (in Southern Italy). This paper reports in Figures 4A,B the kriged probability maps exceeding the erosivity density threshold-value of 3 MJ ha -1 h -1 for the months of August and September, respectively. Our results extend the current knowledge about the hydrological risk prone areas highlighting as ca 65% of the Italian peninsula has a probability higher than 50% to be subjected to hydrological hazards in August. The severely affected hydrological prone areas hazard, in addition to Northern regions already reported by Ballabio et al. (i.e., Friuli-Venezia Giulia, Lombardy, Veneto and Emilia Romagna) are Central regions as Tuscany and Umbria. While Molise, Campania, Basilicata and parts of Sicily are moderately affected (Figure 4A). In the Mediterranean region, circumstances causing precipitation extremes, besides the inherent probability of higher data variance, are due to interaction/combination of different factors acting at local and large scales. Among these, the Mediterranean Sea Surface Temperature (SST), the moisture fluxes from North Atlantic, the orographic forcing factors, and the change of synoptic climatology in different seasons (Figures 4A,B). In particular the characteristic of the development and evolution of cyclones in the Mediterranean have both a smaller spatial and temporal scale than North Atlantic systems, with generally shorter life [START_REF] Trigo | Objective climatology of cyclones in the Mediterranean region[END_REF]. Our results have clearly highlighted the effects on the Italian peninsula, of this strong seasonal signal of the Mediterranean cyclones distribution during precipitation extremes.

In August, Figure 4A reports, over our peninsula, the effect of extreme precipitation events caused by pressure anomalies associated with conditions favoring intense North-Atlantic cyclones (i.e., the cyclogenesis pattern with the cyclone centered North-Western out of the coasts of Ireland and with the low extending south-eastward over the Gulf of Genoa affecting the Alpine region, Tuscany, Umbria regions and the Adriatic sea and the North-Eastern Coast of Sicily. While the extreme precipitation highlighted in September and October (Figures 4B,5A) are associated to the Cyprus low (i.e., the cyclogenesis pattern common in the Mediterranean area during Autumn (SON) and Winter (DJF), which north-westward affects Sicily, Puglia, Ioanian Calabria, Campania and Latium and then the regions northward.

The characteristic depth of Cyprus cyclones, with enhanced (negative) geopotential height anomalies, the intensity ("depth") and its location is the driving factor of the effect of visible heavy precipitations in autumns and winter [START_REF] Maheras | On the vertical structure of composite surface cyclones in the Mediterranean region[END_REF]. The "depth" of a cyclone is the pressure difference between its center and periphery. Our results show this situation. In September, the area hit by storm erosive events (erosivity density >3 MJ ha -1 h -1 ) drops to around 50% of the peninsula. However, the probability of erosive events is still high in along the coastal and insular areas. Veneto, Westerly Umbria, some areas of Liguria, Tuscany and Campania coast still can face hydrological hazard in September (Figure 4B) but the effect of the cyclogenesis characteristics are evident in this month on the coasts and mainland of the Latium, Northern Campania, Puglia Garganica, Ionian Calabria, Sicily and Eastern Sardinia which all are facing higher hydrological hazard as consequence of the Northwestward exstension of the Cyprus Low.

October presents a smaller hydrological hazard, with areas prone to erosive events falling to about 30% of Italian territory; mainly located in the Campania coastal, Southern Basilicata, Sicily, Eastern Sardinia and Latium (Figure 5A).

Taking into account the whole period August-October, the regions Campania, Latium and Sicily are more affected than others. In particular, the months of August and September represent a vulnerable timespan of the year, where both storm are the most erosive and lands are bare. This combination drives to hydrological hazards events across to more than 50% of the Italian territory. Focusing in Latium region, we find that lower Tiber river basin of the Rome metropolitan area (Latina Province) is at a high hazard (Figure 5B). In particular, we can observe the Latina landscape (black square in Figure 5B) concerned in a 3-D view with the erosivity density threshold-value z k of 3 MJ ha -1 h -1 for the October month (Figure 6).

The pink color area in this figure represents the zones where the 3 MJ ha -1 h -1 threshold value has a probability higher than 50%, and thus are the places mostly prone to damaging hydrological hazard at beginning of the autumn. These areas are included among the Villages of Cori, Pontenuovo, Sermoneta, Bassiano, Carpineto Romano and Sezze. This suggests several hydrological implications, which must be taken into account when, for instance, the effects of land use like urbanization and deforestation are examined, or when the relation between climate change and flash-flood occurrence is studied [START_REF] Marchi | Characterisation of selected extreme flash floods in Europe and implications for flood risk management[END_REF].

Vulnerable timespan keeps is very active also for the Sicily island in the month of October. Sicily presents geoclimatic characteristics that make its territory prone to flashflood, especially in September and October, as demonstrated by some recent events [START_REF] Diodato | Storminess and environmental changes in the Mediterranean Central Area[END_REF]Aronica FIGURE 6 | Spatial patterns of kriged-probability perspective view exceeding the erosivity density threshold-value z k of 3 MJ ha -1 h -1 with probability >50% (pink color) for October across the Latina Province (Southern Latium) for the period 2002-2011. Wood areas are in green color. et al., 2012;[START_REF] Fiorillo | Reconstruction of a storm map and new approach in the definition of categories of the extreme rainfall, northeastern Sicily[END_REF][START_REF] Candela | Probabilistic flood hazard mapping using bivariate analysis based on copulas[END_REF]. On the other hand, from the specific analysis of the extreme rainfall occurred in Sicily, it emerges that precipitations are nearly always associated with particularly intense storm [START_REF] Caccamo | Highlights on extreme meteorological events in Sicily[END_REF]. These finding agree with the results of [START_REF] Messeri | Weather-related flood and landslide damage: a risk index for Italian regions[END_REF], who showed the most perturbed flows in Italy, associated with the highest impact in terms of damaging hydrological events, are observed in autumn. These cyclonic perturbations tend to have short life cycles, with average radii in the range 300 to 500 km in this period of the year [START_REF] Lionello | Cyclones in the Mediterranean region: climatology and effects on the environment[END_REF], and a predominance of convective storms. Convective storms can have, however, smaller radii (around 10-50 km), as suggested by the small kriged-ranges obtained with the covariance functions. These are the highest erosive storms that may occur in major part of Italy, characterized by a complex property in transferring erosive energy to lands. In these localities, subgrid scale convection and intensification phenomena dominate rain-producing processes (e.g. [START_REF] Mazzarella | Multifractal dynamic rainfall processes in Italy[END_REF][START_REF] Dünkeloh | Circulation dynamics of Mediterranean precipitation variability 1948-1998[END_REF], accompanied by several rain episodes releasing in a few hours as much energy as in an average year or more. These changes capture the shift toward autumn of some flash-flood conditions typical of summer [START_REF] Gaume | A compilation of data on European flash floods[END_REF]. Delayed effects show that the warmer Mediterranean period taking place from end-summer to early autumn is leading to more frequent and intense torrential rains [START_REF] Millán | Climatic feedbacks and desertification: the Mediterranean model[END_REF].

Cross-Validation and Error Assessment

The results from the cross-validation and error assessment procedures check have been presented with joint statistical errors and as monthly scatter diagrams (Figures 7A-C). They show that the measured values of erosivity density are in agreement with the predicted probability indicators.

In this way, mean values of 0.006 in August and September, and 0.01 in October, make evident that there is no systematic error because the root mean square error is close to the average standard errors. The mean thus correctly assesses the variability in predictions. In addition, the Root Mean Square Standardized (RMSS) error variance is the same and close the 1.0 as compared to the theoretical variance of the kriging variance.

At sub-regional scale, we performed a qualitative error assessment by comparing kriged-probability erosivity density >3 MJ ha -1 h -1 (Figure 8A) to previously published damaging hydrological risk maps (joint probability of landslides and floods, Figure 8B, [START_REF] Guzzetti | An Information System on Historical Landslides and Floods in Italy. Urban Hazards Forum[END_REF] for the Sicily island.

The comparison of kriged probability erosivity density higher than the threshold (orange colors in Figure 8A) is in agreement with probability hydrological risk areas (blue colors in Figure 8B), except for the extreme westernmost part of the island. However, a strict parallelism between the two maps is not to be expected considering the multifaceted nature of the involved phenomena and the different time periods covered.

CONCLUSION

The widespread availability of high-temporal resolution rainfall records for large areas and the development of climate models opens new opportunities for using geostatistical methods for large scale planning and hazard prevention. This paper presents a geostatistical modeling framework, which uses a probabilistic approach to assess the spatial variability of storm erosivity density throughout the Italian regions. The results so obtained suggest the possibility of using geostatistical spatial modeling to determine the probability of exceeding high erosivity density thresholds and generate probability maps to delineate the most sensitive areas. These areas, highly vulnerable, may result in catastrophic regime shifts connected with the occurrence of damaging hydrological events. We offer the results of this work as a springboard to support policy-makers, ) for October with (B) damaging hydrological risk maps (joint probability of landslides and floods occurrence) with 5-years return period [START_REF] Guzzetti | An Information System on Historical Landslides and Floods in Italy. Urban Hazards Forum[END_REF] for the Sicily region.

local authorities and civil protection in planning actions, in the medium and long term, aimed at reducing hydrological disasters. In fact, such a soft-computing model represents a paradigm shift on how to provide timely, accurate and actionable hydrologic hazard information. Relying on credible information regarding erosivity density and its environmental drivers, this study shows how geostatistical methods can be practically implemented to create countrywide spatially explicit probabilistic maps capable to describe the seasonal effect of cyclogenesis' modification as the case of the highlighted seasonal variation of the North Atlantic and Cyprus Low systems during extreme events. These probabilistic maps can therefore be of great help in studying erosive hazards or the effects over Italy, of the generation and growth of Mediterranean cyclones, which however depends on the relationship of several factors among others, the complex orography, large scale air flow, and SST. In addition, these probabilistic maps refer not only to a possible check of the physical understanding of some processes, but also to the possibility to gain information where data or measurements actually do not exist. In fact, without giving the value in each point but returning a map of probabilities, this approach can support decision-making because it identifies the hydrologic danger associated with the probability to exceed an erosivity density threshold.

ED-outputs have the advantage to be developed based on more available data and can be scaled up for producing weather hazard maps. Vice versa, the disadvantage is that the same probability to overcome a specific ED threshold may have different meanings when superimposed to areas with different combinations of morphology, soil erodibility, land use, plant density, and management. In addition, these findings are purely geostatistical, with no explicit mechanistic modeling of erosivity density.

This study has also some limitations because it does not consider that population density, infrastructures, plant density, and other factors also influence the occurrence of a damage. In addition, the particularly prolonged low-intensity rainy periods, which can condition landslide fatalities, are not accounted in this probability calculation.

FIGURE 1 |

 1 FIGURE 1 | (A) Overview map of the study area. (B) The spatial pattern of mean annual rainfall amount overlaid with stations location, and (C) Ratio of maximum-recorded daily rainfall vs. mean annual precipitation for Italy during the period 2002-2011.

  has been used. Stage 1 assumes an isotropic and spherical model and computes empirical covariance and cross-covariance functions on standard deviationscaled data. A covariance distance measures the average degree of dissimilarity between unsampled values and nearby values. Empirical indicator covariance function C X (h; z k ) and cross-covariances CC XI (h; z k ) are thus produced (see section Results and discussion). Stage 2, calibrated interactively the geostastistical parameters wich include: range (a), number of lag (assumed equal 7), or lag size h (assumed equal 24 km), model type, nugget and partial sill.

FIGURE 2 |

 2 FIGURE 2 | Isotropic neighboring circle for erosivity density during (A) August, (B) September, and (C) October. The colured bands in the corner represent the patterns of station-points inside of any neighboring circle with the related percentage weighted given for any data-points.

FIGURE 3 |

 3 FIGURE 3 | Covariance function C X (h) at threshold z k (ED>3 MJ ha -1 h -1 ) for (A) August, (B) September, and (C) October month, with respective model of co-regionalization between indicator and uniform transforms [cross-covariance CC XI (h; z k )], (a 1 ), (b 1 ), and (c 1 ).

FIGURE 4 |

 4 FIGURE 4 | Spatial patterns of kriged-probability maps exceeding the erosivity density threshold-value z k of 3 MJ ha -1 h -1 for (A) August and (B) September across Italy for the period 2002-2011.

FIGURE 5 |

 5 FIGURE 5 | Spatial patterns of kriged-probability maps exceeding the erosivity density threshold-value z k of 3 MJ ha -1 h -1 for (A), October across Italy for the period 2002-2011 and (B) zoomed area in Central Italy where boundaries of river basin and main cities are over-imposed. The rectangular shape indicates the area for which the perspective view of Figure 6 was made.

FIGURE 7 |

 7 FIGURE 7 | Scatterplot of cross-validation between storm erosivity measured and indicator probability kriging for the erosivity-thresholds z k (ED>3 MJ ha -1 h -1 ), for (A) August, (B) September, and (C) October month. The blue vertical lines in both the graphs represent the respective ED thresholds, while the gray horizontal lines separate the kriged-prediction indicator with probability above 0.50 from those below 0.50.

FIGURE 8 |

 8 FIGURE 8 | Comparison of (A) kriged probability for erosivity density (ED>3 MJ ha -1 h -1 ) for October with (B) damaging hydrological risk maps (joint probability of landslides and floods occurrence) with 5-years return period[START_REF] Guzzetti | An Information System on Historical Landslides and Floods in Italy. Urban Hazards Forum[END_REF] for the Sicily region.

Frontiers in Environmental Science | www.frontiersin.org 2 December 2019 | Volume 7 | Article 193

Frontiers in Environmental Science | www.frontiersin.org 3 December 2019 | Volume 7 | Article 193

Frontiers in Environmental Science | www.frontiersin.org

December 2019 | Volume 7 | Article 193

Frontiers in Environmental Science | www.frontiersin.org 6 December 2019 | Volume 7 | Article 193

Frontiers in Environmental Science | www.frontiersin.org 9 December 2019 | Volume 7 | Article 193

Frontiers in Environmental Science | www.frontiersin.org 10 December 2019 | Volume 7 | Article 193

ACKNOWLEDGMENTS

We thank the NTNU publishing fund support to cover article processing charges in Open Access Journals.

DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS

ND conceived the study, performed the analysis and drafted the manuscript with GB, who wrote the final manuscript. PB and PP prepared and organized the sediment data used in the study, and contributed to scientific discussion of the article together with CB.

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.