# QPF performance of the updated SIMM forecasting system using reforecasts

Stefano Mariani, Marco Casaioli, Alexandre Lanciani, Simona Flavoni, Christophe Accadia

<figure><img src="image_1.png" /></figure>

To cite this version:

Stefano Mariani, Marco Casaioli, Alexandre Lanciani, Simona Flavoni, Christophe Accadia. QPF performance of the updated SIMM forecasting system using reforecasts. Meteorological Applications, 2015, 22 (2), pp.256 - 272. 10.1002/met.1453 . hal-01496783

HAL Id: hal-01496783

https://hal.science/hal-01496783v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# QPF performance of the updated SIMM forecasting system using reforecasts

Stefano Mariani,<sup>a,*</sup> Marco Casaioli,<sup>a</sup> Alexandre Lanciani,<sup>b</sup> Simona Flavoni<sup>c</sup> and Christophe Accadia<sup>d</sup>

**ABSTRACT:** In late 2009, the meteorological component of the Institute for Environmental Protection and Research (ISPRA)'s Hydro-Meteo-Marine Forecasting System (*Sistema Idro-Meteo-Mare* – SIMM), namely the hydrostatic Quadrics BOlogna Limited Area Model (QBOLAM), was replaced by an up-to-date version of the same model. This new version, which is simply referred to as BOLAM, deploys advanced parameterization schemes for cumulus convection, radiation, soil, and exchange of turbulent fluxes and uses the explicit advection of five hydrometeors. As a part of the SIMM verification programme, it was judged objectively how these changes affected the quality of quantitative precipitation forecasts (QPFs). This was done by intercomparing the current operational BOLAM against QBOLAM. Observational analyses, from April to September 2001, deployed earlier to evaluate the QBOLAM performance, were considered for this intercomparison. An *ad hoc* reforecast campaign was performed to generate the database of the corresponding retrospective QPFs for BOLAM. A multimethod verification approach was used to assess the spatial scales resolved by the intercompared model versions and the skill at forecasting both frequent, high-base rate and rare, low-base rate precipitation events. The effect on performance measures of differences in bias between the intercompared versions was also evaluated. The results indicate an improved performance of BOLAM compared with the older version not only when considering the native 0.1° model grid, where the intercomparison may be affected by differences in the resolved scales, but also when using a coarser 0.5° grid, where the two remapped QPF series show similar spatial scales. Fewer false alarms are produced by BOLAM during the verification time period.

**KEY WORDS** forecast verification; intercomparison; power spectrum analysis; categorical scores; high-base rate and low-base rate precipitation events; sensitivity to bias of performance measures; hydrostatic limited area model

Received 10 May 2013; Revised 18 December 2013; Accepted 6 January 2014

## 1. Introduction

The improvement of weather forecasts is one of the primary goals of any meteorological, hydrological or environmental institution running a numerical weather prediction (NWP) model. Weather predictions can be improved by enhancing the NWP model accuracy, for instance by implementing more accurate and advanced physical parameterizations, improving the quality of the initial and boundary conditions (IC/BCs), or increasing the horizontal and vertical resolutions (see, e.g. Brooks and Doswell, 1993; Gallus and Segal, 2001; Mass et al., 2002; Hamill et al., 2004, and references therein). Porting the NWP model onto a more efficient computer machine might also have an effect on its forecasting performance. To evaluate statistically in what way and how much such changes affect the model performance, it becomes essential to verify (over a sufficiently long period) the weather forecasts generated by the updated model version and, in addition, to compare the verification results with those obtained from previous model

versions. The focus of the present work is the assessment of model changes on quantitative precipitation forecasts (QPFs).

The availability of computational power capabilities has made it possible to generate extended datasets of retrospective forecasts (reforecasts) using a NWP model in relatively short times. More generally, reforecast datasets have the appeal to be useful either to diagnose model biases, to correct statistically weather forecasts, or to study the atmospheric predictability (e.g. Hamill et al., 2006, 2013; Hagedorn et al., 2012). For diagnosing improvements in an updated forecasting system, it can be a natural choice for the model developers to generate the retrospective forecasts over the same period(s) used for verifying previous system versions in order to prevent different datasets of observations (against which forecasts are evaluated) from affecting the comparison of the verification results.

The choice of the verification methodology used for inter-comparing the different model versions is also important. Usually, verification is based on performing exploratory data analysis and computing one or more performance measures (referred to as traditional or standard verification methods – see WWRP/WGNE Joint Working Group on Verification, 2004; Wilks, 2006; Jolliffe and Stephenson, 2011; WWRP/WGNE Joint Working Group on Forecast Verification Research, 2013, and references therein), which evaluate the correspondence between forecasts and observations. In the last two decades,

<footnote><sup>a</sup> Institute for Environmental Protection and Research (ISPRA), Water Protection Department, Rome, Italy</footnote>

<footnote><sup>b</sup> IFS Italia, Rome, Italy</footnote>

<footnote><sup>c</sup> Laboratoire d'Océanographie et du Climat: Expérimentation et Approches Numériques, Institut Pierre Simon Laplace (LOCEAN-IPSL), Unité Mixte de Recherche 7617 CNRS, Paris, France</footnote>

<footnote><sup>d</sup> European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT), Remote Sensing and Products Division, Darmstadt, Germany</footnote>

<footnote>* Correspondence: S. Mariani, Institute for Environmental Protection and Research (ISPRA), Via Vitaliano Brancati 48, 00144 Rome, Italy. E-mail: stefano.mariani@isprambiente.it</footnote>
though, spatial verification methods have become more and more popular in the forecast verification community, because they provide thorough information on forecast quality (see Davis et al., 2006a, 2006b; Casati et al., 2008; Ebert, 2008; Gilleland et al., 2009, 2010, and references therein). On one hand, weather forecasters and model developers appreciate the ability of these new methods to mimic computationally the human judgement in assessing the forecast accuracy and to quantify and qualify the sources of forecast error (see, e.g. Hoffman et al., 1995; Ebert and McBride, 2000; Grams et al., 2006; Moise and Delage, 2011). On the other hand, as pointed out by Rossa et al. (2008), maintaining the calculation of traditional measures is worthwhile when verification results related to previous forecasting systems (old model versions) are available and they can be used as benchmarks to assess the accuracy of a newer forecasting system.

This was the premise of an extensive verification programme, undertaken between mid-2007 and mid-2012, aimed at evaluating the quality of the forecasts generated by the meteorological component of the ISPRA's Hydro-Meto-Marine Forecasting System (Sistema Idro-Meto-Mare - SIMM; Speranza et al., 2004, 2007; Mariani et al., 2005). Over the past 6 years, several forecast parameters, including geopotential height, mean sea level pressure, potential vorticity, precipitation, total column water vapour, and surface wind, were evaluated over single intense/severe case studies and over long time periods. In either case, traditional verification measures and feature-based methods were deployed in verifying the SIMM meteorological forecasts. One of the main aims was to assess the impact on QPFs of the physical parameterization schemes introduced into an up-to-date parallel version of the hydrostatic BOlogna Limited Area Model (BOLAM - Buzzi et al., 1994; Malguzzi and Tartaglione, 1999), developed by the Institute for Atmospheric Sciences and Climate of the Italian National Research Council (ISAC-CNR), that in late 2009 replaced within SIMM, a previous version named Quadrics BOLAM (QBOLAM). In order to make the QPF verification results of this newer BOLAM version fully comparable to those related to the former implemented model versions, a BOLAM reforecast dataset was generated over a 6 month period, April to September 2001, exploited for previous studies (Accadia et al., 2005; Mariani, 2008). Hence, the same precipitation objective analyses used in those studies, which were obtained from rain gauge measurements covering the entirety of Italy, were considered for the comparison. The choice was also connected to the fact that QBOLAM exhibited the worst performance during this time period, mostly owing to a systematic overprediction of the occurrence of intense precipitation events.

Following the rationale proposed by Rossa et al. (2008) the QPF verification was based mainly on the calculation of several traditional categorical scores, the same ones deployed in the aforementioned studies. In addition, newly proposed scores for evaluating intense precipitation forecasts and assessing the sensitivity to bias of the competing model versions were used. Bootstrap-based confidence intervals were calculated to assess statistically the differences in the scores. Spatial methods were not applied in this work. The application of such methods, which are usually much more suitable for higher resolution verification studies, may be problematic when dealing with small or narrow verification areas, as in the case of Italy (see, e.g. Tartaglione et al., 2008). The spatial distributions of the contingency table frequencies used to calculate the categorical scores were, however, investigated to assess whether the forecast performance of the different

model versions showed location dependencies. The forecast verification was also achieved through the use of the power spectrum analysis of QPFs to relate the differences in scores to the differences in the spatial scales resolved by the compared versions.

This paper is organized as follows. Section 2 provides a brief overview of the structure of the SIMM forecasting system. QPF fields and observational rainfall analyses are explained in Section 3. Section 4 describes the multimethod approach used in this intercomparison study. Results are discussed in Section 5. Finally, Section 6 provides conclusions and final remarks on future activities.

## 2. The SIMM forecasting system

### 2.1. The modelling chain

SIMM is a meteo-marine forecasting system that is operational at ISPRA (formerly APAT) since 2000. This modelling chain was primarily designed in the late 1990s to provide an integrated system over the entire Mediterranean area, which is able to combine analysis and forecast of weather together with the relevant sea surface and hydrological post-processes (Speranza et al., 2007). The main core of SIMM is formed by a cascade of four numerical models (Figure 1) telescoping from the Mediterranean basin to the Venice Lagoon. IC/BCs are provided on a 0.5° grid by analyses and forecasts from the 1200 UTC run of the European Centre for Medium-Range Weather Forecasts (ECMWF) global model. To allow a sufficiently smooth transition from the coarser resolution of the 0.5° ECMWF data down to the 0.1° model grid covering the Mediterranean basin (referred to as very high-resolution grid - VHR grid; Figure 2), the NWP model at the base of the chain is one-way nested in two steps using a grid of intermediate extension with horizontal spacing equal to 0.3° (referred to as high-resolution grid - HR grid; Figure 2). For technical reasons, this model design has been kept unchanged until now. ECMWF IC/BCs are now operationally available

<figure><img src="image_4.png" /><figcaption>Figure 1. The SIMM model sequence and operational chain: from the initialization with the ECMWF fields as IC/BCs to the storage and operational use of the weather, marine and coastal forecast products (see Section 2.1 for model acronyms).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 2. Computational domains for the 0.3° HR BOLAM grid (long-dashed line) and the 0.1° VHR BOLAM grid (solid line), with spatial distribution of the 1537 rain gauges (grey points) available over the entire Italy.</figcaption></figure>

for ISPRA on a 0.25° grid, so that an improved SIMM design, which was tested recently (Casaioli *et al.*, 2013), is scheduled to be operational by the beginning of 2014.

The 0.1° BOLAM is then coupled with: (1) a 0.1° WAve Model (WAM – The Wamdi Group, 1988) to simulate the sea state over the Mediterranean Sea; (2) a high-resolution shallow-water version of the Princeton Ocean Model (POM – Blumberg and Mellor, 1987) to forecast the sea surface elevation of the Adriatic and Ionian Sea; and (3) a finite-element hydrodynamic shallow-water model (VL-FEM – Umgiesser and Bergamasco, 1993) aimed at forecasting the exceptional tide peak events (*acqua alta*) that take place periodically in the Northern Adriatic Sea and, in particular, in the Lagoon of Venice.

The SIMM NWP model is also coupled, but only in a research configuration, with the physically based, fully distributed, rainfall-runoff TOPographic Kinematic AProxi-mation and Integration (TOPKAPI – Todini and Ciarapica, 2001). Furthermore, more recently, two new operational com-ponents have been coupled to BOLAM (see Figure 1), namely a Mediterranean Coastal Wave Forecasting System (MC_WAF – Inghilesi *et al.*, 2012) to simulate waves in sev-eral Italian coastal areas and a new storm surge forecasting system, based on the Shallow water HYdrodynamic Finite Ele-ment Model (SHYFEM – Umgiesser *et al.*, 2004) developed at the Institute of Marine Sciences of CNR, to model the sea sur-face elevation in the Mediterranean Sea and the Venice Lagoon.

## 2.2. The evolution within SIMM of the BOLAM model

The NWP model that initially underpinned the SIMM forecast-ing system was QBOLAM, a model version parallelized on the Alenia-Space Co. Quadrics supercomputer, which was based on the BOLAM code dated 1997. Because of the structure of Quadrics, some simplified physical parameterization schemes were adopted, including the Kuo scheme (Kuo, 1974; Raymond and Emanuel, 1993) for cumulus convection. This version, hereinafter referred to as QBOLAM-Kuo, showed an overall fair performance when compared with other operational hydro-static and non-hydrostatic NWP models (Accadia *et al.*, 2003a, 2005; Mariani *et al.*, 2005, 2009). Nevertheless, a systematic tendency was observed for this model version to overpredict intense precipitations during the summertime, especially over

mountain ranges. The use of a simplified cumulus convection parameterization scheme was hypothesized as a possible cause of this shortcoming.

In 2007, the porting of the SIMM forecasting system onto the Silicon Graphics Inc. Altix 8-processor machine provided the opportunity to validate this hypothesis. It was possible to implement into QBOLAM, even if in a research configura-tion, the Kain–Fritsch cumulus convection scheme (Kain and Fritsch, 1990, 1993). This research version, which is hereinafter defined as QBOLAM-KF, was verified against the same series of observational precipitation analyses used for the assessment of the QPF performance of QBOLAM-Kuo (Mariani, 2008). The results showed a statistically significant QPF improve-ment, with a sensible reduction of false alarms, especially over the summer season, when moving from QBOLAM-Kuo to QBOLAM-KF, confirming a lesser capability of the Kuo parameterization scheme to produce reliable high-resolution (0.1°) rainfall fields.

Following this outcome, it was decided not only to implement operationally the Kain–Fritsch parameterization scheme but also to investigate further the use of an updated BOLAM version into the SIMM chain. A fully up-to-date BOLAM version (identified here as BOLAM-KF), which is the core of this statistical analysis of comparative performance, became operational into SIMM in late 2009. BOLAM-KF deploys the 3D Eulerian-weighted averaged flux scheme for advection, the explicit advection of five hydrometeors (cloud ice, cloud water, rain, snow, and hail/graupel), the Schultz microphysical scheme for large-scale precipitation, the Kain–Fritsch parameterization scheme for cumulus convec-tion, a combination of the Geleyn and Morcrette schemes for radiation, the E-l scheme for exchange of turbulent fluxes, and a four-layer advanced scheme for soil (for more details, see http://www.isac.cnr.it/dinamica/projects/forecasts/bolam_short_description_2012.htm).

## 3. QPF fields and observational rainfall analyses

Two precipitation reforecast series for the period April to September 2001 (183 days) were considered for the intercom-parison: one obtained by joining together the daily retrospective
<table><caption>Table 1. Schematic description of the principal model settings of QBOLAM-KF and BOLAM-KF.</caption><thead><tr><th colspan="6">Initialization settings</th></tr><tr><th>HR model version</th><th>VHR model version</th><th>HR model IC</th><th>HR model BC</th><th>VHR model IC</th><th>VHR model BC</th></tr></thead><tbody><tr><td>0.3° QBOLAM-KF</td><td>0.1° QBOLAM-KF</td><td>1200 UTC ECMWF analysis</td><td>1200 UTC ECMWF forecast every 6 h</td><td>0.3° QBOLAM-KF nesting</td><td>0.3° QBOLAM-KF every 3 h</td></tr><tr><td>0.3° BOLAM-KF</td><td>0.1° BOLAM-KF</td><td>1200 UTC ECMWF analysis</td><td>1200 UTC ECMWF forecast every 6 h</td><td>0.3° BOLAM-KF nesting</td><td>0.3° BOLAM-KF every 3 h</td></tr></tbody><thead><tr><th colspan="6">Timing settings</th></tr><tr><th>Model version</th><th>HR model initial time</th><th>Spin-up time</th><th>VHR model initial time</th><th>VHR model forecast range<sup>a</sup></th><th>Output frequency</th></tr></thead><tbody><tr><td>QBOLAM-KF</td><td>1200 UTC of the previous day</td><td>12 h</td><td>0000 UTC</td><td>24 h</td><td>3 h</td></tr><tr><td>BOLAM-KF</td><td>1200 UTC of the previous day</td><td>12 h</td><td>0000 UTC</td><td>48 h</td><td>1 h</td></tr></tbody><thead><tr><th colspan="6">Forecasting settings</th></tr><tr><th>Model version</th><th>Advection</th><th>Cumulus convection</th><th>Radiation</th><th>Turbulence</th><th>Soil</th></tr></thead><tbody><tr><td>QBOLAM-KF</td><td>Forward-backward advection scheme</td><td>Kain-Fritsch scheme</td><td>Page scheme</td><td>Louis scheme</td><td>2 levels + 1 (climatological) level</td></tr><tr><td>BOLAM-KF</td><td>Weighted average flux scheme</td><td>Kain-Fritsch scheme</td><td>Geleyn and Morcrette schemes (staggered)</td><td>E-l scheme</td><td>3 levels + 1 (climatological) level</td></tr></tbody></table>

forecasts generated with the current operational 0.1° QBOLAM-KF version and the other obtained by joining together the daily retrospective forecasts generated with the 0.1° QBOLAM-KF version (the analysis benchmark). A schematic description of the settings used to generate these two reforecast series is reported in Table 1.

For the study purposes, the daily reforecasts were obtained by accumulating the precipitation fields modelled from 0000 UTC to 2400 UTC (+12 h to +36 h from the ECMWF initialization time). Intercomparison was limited here to this short lead time, because reforecasts for the 0.1° QBOLAM-KF were produced only for 24 h (see Table 1). Given that the current operational forecast range for the 0.1° QBOLAM-KF is 84 h (see http://www.isprambiente.gov.it/pre_meteo/simm.html; in Italian), a further verification study is planned to evaluate the model performance for longer lead times.

One of the controversial questions that arises when verifying forecasts concerns the approach that should be used to treat observations. Several researchers in recent years have faced this question, including Skelly and Henderson-Sellers (1996), Ghelli and Lalaurette (2000), Ebert (2001), Cherubini et al. (2002), Ebert et al. (2003), WWRP/WGNE Joint Working Group on Verification (2004), Ahrens and Beck (2008), Göber et al. (2008) and Jolliffe and Stephenson (2011). The so-called grid point approach, which consists of interpolating forecasts onto the station locations, is usually opposed to the so-called grid box approach, which consists of up-scaling the observations by computing an observational analysis over the model grid or, more generally, over a coarser verification grid. Model type, physical variable analysed, and forecast verification purpose can affect the choice of using one approach rather than the other. The grid box approach can also be considered more model-oriented, whilst the comparison at gauges can be valid, and necessary, from the user's perspective (Ebert, 2001).

Nonetheless, both present pros and cons, including the introduction of a 'representativeness' error (Tustison et al., 2001) due to representing either forecasts (in the grid point approach) or observations (in the grid box approach) to a different scale

compared with the own inherent one. The effectiveness of adopting a grid box approach when verifying QPF has been proved by Ghelli and Lalaurette (2000), Cherubini et al. (2002), Ahrens and Beck (2008), and Göber et al. (2008), even when observations are available over alpine or mixed area, as in the case of Italy. On the basis of these outcomes, observational rainfall analyses are typically deployed in ISPRA to evaluate the grid box behaviour described by the SIMM meteorological model. This choice was also retained in the present work to maintain the same observational reference of previous studies, thus making comparable the results obtained with the different model versions in this manner.

The daily rain gauge analyses used in the intercomparison were obtained as follows. First, precipitation measurements were collected from the former network of the Italian National Hydrographic and Marigraphic Service (SIMN, now at ISPRA) and from the regional networks of Calabria, Emilia-Romagna, Liguria, Marche, Piedmont, Sardinia, Sicily and Valle d'Aosta regions. A total of 1183 of the 1537 rain gauges available (Figures 2 and 3) were working during the considered period and around two-thirds of these gauges had an 80-100% temporal coverage, guaranteeing a good sample dataset of measurements across Italy. These rainfall measurements were then quality checked, accumulated on a daily basis from 0000 UTC to 2400 UTC and, finally, the so-obtained 24 h precipitation values were analysed objectively using a two-pass Barnes scheme (Barnes, 1964, 1973; Koch et al., 1983; Accadia et al., 2005; Lanciani and Salvati, 2008). This is a Gaussian-weighted averaging interpolation scheme whose weights are a function of the averaged spacing among the considered gauges and the distance between a gauge and the centre of the verification grid box. The first pass of the Barnes scheme produces a first-guess precipitation analysis, which is then followed by the second pass to increase the amount of details from the prior pass.

In this study, the convergence parameter was equal to 0.3, whereas the average data spacing was set to 0.2°. This latter setting fulfilled the constraint that the ratio between the model grid size and the average data spacing lies in the range of

<footnote><sup>a</sup>The forecast range indicated refers to the one used to generate the retrospective forecast time series.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 3. Spatial distribution of the 1183 rain gauge stations working during April to September 2001 plotted as a function of the 99th percentile precipitation value: blue points indicate values ≤ 20.0 mm 24 h⁻¹; red points indicate values between 20.0 and 40.0 mm 24 h⁻¹; green points indicate values > 40.0 mm 24 h⁻¹.</figcaption></figure>

0.3–0.5 (Barnes, 1964, 1973). Grid points that did not have any rain gauge within a radius of 0.15° were neglected from the analysis in order to avoid the excessive rainfall spreading introduced by the Barnes scheme on model grid points far from the actual locations of rain gauges. A rough estimate of the uncertainty associated with this scheme was obtained by calculating the difference between the Barnes gridded analysis and a gridded field obtained by averaging all the observations within each grid box. The difference averaged over the whole verification time period ranges from –6.4 to 3.6 mm, displaying the maximum discrepancy in very few localized areas characterized by intense precipitations (not shown). On average, the discrepancies are of the order of 3.0 mm and they tend to be concentrated in areas where strong horizontal gradients of rainfall amount are present. These differences are considered as negligible for the purposes of this statistical verification study.

# 4. The verification multimethod approach

## 4.1. Categorical scores

The calculation of the categorical scores is based on a binary ‘yes/no’ statement by considering whether the QPFs and the precipitation objective analyses exceed or not a range of predefined thresholds. In the present study, these rainfall thresholds are set to 0.5, 5.0, 10.0, 20.0, 30.0 and 40.0 mm 24 h⁻¹. This choice is consistent with previous settings (see Accadia *et al.*, 2005; Mariani, 2008) that took into consideration the observed precipitation distribution for the definition of the high-base rate events, and the distribution of the 99th percentile precipitation amounts for the definition of the rare, low-base rate events (see Figure 3). The possible combinations of these dichotomous statements are summarized in the 2 × 2 contingency table (see Table 2), whose elements represent the

<table><caption>Table 2. The 2 × 2 contingency table, representing the combination of the event possibilities of a sequence of N observation–forecast pairs, expressed in the top in terms of hits (a), false alarms (b), misses (c), and correct non-rain forecast (d) and in the bottom in terms of quantity of hits (H), forecasts of the event (F), observations of the event (O), and N.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">Rain event observed</th></tr><tr><th>Yes</th><th>No</th><th>Total</th></tr></thead><tbody><tr><th rowspan="3">Rain event forecast</th><th>Yes</th><td>a</td><td>b</td><td>a + b</td></tr><tr><th>No</th><td>c</td><td>d</td><td>c + d</td></tr><tr><th>Total</th><td>a + c</td><td>b + d</td><td>a + b + c + d = N</td></tr></tbody></table><table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">Rain event observed</th></tr><tr><th>Yes</th><th>No</th><th>Total</th></tr></thead><tbody><tr><th rowspan="3">Rain event forecast</th><th>Yes</th><td>H</td><td>F - H</td><td>F</td></tr><tr><th>No</th><td>O - H</td><td>N - F - O + H</td><td>N - F</td></tr><tr><th>Total</th><td>O</td><td>N - O</td><td>N</td></tr></tbody></table>

observed rain events detected correctly (a or ‘hits’), the rain events forecasted but not observed (b or ‘false alarms’), the observed rain events not detected (c or ‘misses’), and the events neither observed nor forecasted (d or ‘correct non-rain forecasts’). The contingency table frequencies can also be expressed in terms of quantity of hits (H = a), forecasts of the event (F = a + b), observations of the event (O = a + c), and total number of the compared observation–forecast pairs (N = a + b + c + d).

Several scores have been developed, introduced and used in the forecast verification to evaluate the NWP model performance accurately. However, different scores assess different attributes of the forecast quality (Murphy, 1993); hence, in this study, as suggested by Stephenson (2000) and Nurmi (2003), it was decided to apply more than one verification measure. The scores considered (described in Table 3) are the frequency bias (BIAS), the probability of detection (POD), the false alarm ratio (FAR), the equitable threat score (ETS), the Hanssen–Kuipers score (HK), the extreme dependency index (EDI), and the symmetric extreme dependency index (SEDI). The two latter indices are introduced to verify the forecast accuracy in correctly predicting the events observed with a low probability, because they have the appealing propriety of not degenerating to trivial values at high precipitation thresholds.

A forecast system is unbiased when the BIAS score is equal to 1, whereas a value greater (less) than 1 indicates that the forecast system tends to overpredict (underpredict) the occurrence of precipitation events. For a perfect forecast system, ETS, HK, POD, EDI and SEDI are equal to 1 and FAR is equal to 0. A forecast system with no skill shows values equal to 0 or negative (depending on the score range; see Table 3) for ETS, HK, POD, EDI and SEDI and a value equal to 1 for FAR.

As emphasized by Ferro and Stephenson (2011), the calculation of EDI and SEDI requires recalibrating the investigated QPF series to remove the effect of bias. This was done by introducing in the 2 × 2 contingency table different, adjusted thresholds for forecasts so as to guarantee recalculated BIAS values equal to 1 up to a small irreducible absolute error term, εb, arbitrarily set here to 0.05. Computationally, these thresholds were selected by applying the BIAS adjustment method introduced by Hamill (1999) to each one of the two model
<table><caption>Table 3. Summary of verification measures deployed in the study.</caption><thead><tr><th>Name of measure</th><th>Definition</th><th>Range</th></tr></thead><tbody><tr><td>Frequency bias (BIAS)</td><td>BIAS = <sup>a+b</sup>&frasl;<sub>a+c</sub></td><td>[0,∞]</td></tr><tr><td>Equitable treat score (ETS)</td><td>ETS = <sup>a-a<sub>t</sub></sup>&frasl;<sub>a+b+c-a<sub>t</sub></sub> = <sup>H-FO/N</sup>&frasl;<sub>F+O-H-F</sub>, where a<sub>t</sub> = (<sup>a+b</sup>&frasl;<sub>a+c</sub>)(a+c) = <sup>FO</sup>&frasl;<sub>N</sub> are the hits expected from a random forecast</td><td>[-1/3, 1]</td></tr><tr><td>Hanssen-Kuipers score (HK)</td><td>HK = <sup>ad-bc</sup>&frasl;<sub>(a+c)(b+d)</sub></td><td>[-1, 1]</td></tr><tr><td>False alarm ratio (FAR)</td><td>FAR = <sup>b</sup>&frasl;<sub>a+b</sub></td><td>[0,1]</td></tr><tr><td>Probability of detection (POD)</td><td>POD = <sup>a</sup>&frasl;<sub>a+c</sub></td><td>[0,1]</td></tr><tr><td>Probability of false detection (POFD)</td><td>POFD = <sup>b</sup>&frasl;<sub>b+d</sub></td><td>[0,1]</td></tr><tr><td>Bias-adjusted ETS using dH/dF (ETSA<sub>dH/dF</sub>)</td><td>ETSA<sub>dH/dF</sub> = <sup>H<sub>a</sub>-O<sup>2</sup>/N</sup>&frasl;<sub>O(2-O/N)-H<sub>a</sub></sub>, where F = O and H = H<sub>a</sub> = O(1 - (<sup>O-H</sup>&frasl;<sub>O</sub>)) is obtained via dH/dF ∝ (O - H)</td><td>[-1/3, 1]</td></tr><tr><td>Bias-adjusted ETS using dH/dA (ETSA<sub>dH/dA</sub>)</td><td>ETSA<sub>dH/dA</sub> = <sup>H<sub>a</sub>-O<sup>2</sup>/N</sup>&frasl;<sub>O(2-O/N)-H<sub>a</sub></sub>, where F = O and H = H<sub>a</sub> = O - (<sup>O-H</sup>&frasl;<sub>O-H</sub>) Lambert W (<sup>O</sup>&frasl;<sub>O-H</sub> ln(<sup>O</sup>&frasl;<sub>O-H</sub>)) is obtained via dH/dA ∝ (O - H) with dA = dF - dH and Lambert W(.) the omega/Lambert W function</td><td>[-1/3, 1]</td></tr><tr><td>Extreme dependency index (EDI)</td><td>EDI = <sup>log(POFD)-log(POD)</sup>&frasl;<sub>log(POFD)+log(POD)</sub></td><td>[-1, 1]</td></tr><tr><td>Symmetric extreme dependency index (SEDI)</td><td>SEDI = <sup>log(POFD)-log(POD)-log(1-POFD)+log(1-POD)</sup>&frasl;<sub>log(POFD)+log(POD)+log(1-POFD)+log(1-POD)</sub></td><td>[-1, 1]</td></tr></tbody></table>

versions under investigation after considering a hypothetical unbiased forecast system as 'reference'. This method consists of introducing adjusted forecast thresholds for the 'competitor' forecast system that are selected appropriately so as to obtain adjusted BIAS values close to those of the 'reference' forecast system. Another reason for applying the Hamill approach in this study was to remove the bias dependency from the result of the score comparison. Erroneous or ambiguous assessment of the differences in performance scores between competing forecasts may be the consequence of relative high differences in the frequency bias values. It is not unusual that, under particular conditions, a biased forecast system (despite its large number of false alarms) scores better than an unbiased (or a less biased) forecast system, because of the tendency of the former to over-predict the occurrence of rain events (see, e.g. Shuman, 1980; Mason, 1989; Hamill, 1999). Hence, the ETS, HK, POD and FAR scores were also recalculated using, for the more biased forecast system (the 'competitor'), adjusted forecast thresholds achieved in a manner to ensure BIAS values consistent with the 'reference' BIAS values.

In addition, the sensitivity to bias of performance measures was accounted for by means of the two bias-adjusted ETSs introduced by Mesinger and Brill (2004) and Mesinger (2008). These scores, named ETSA<sub>dH/dF</sub> and ETSA<sub>dH/dA</sub> (described in Table 3), correct the equitable threat for the bias dependency by extrapolating or interpolating H to the value that the forecasting system would have if it had no bias (BIAS = F/O = 1). These scores, filtered out by the dependence on bias, provide a measure of the forecast accuracy in placing the precipitation above a considered threshold. This bias-adjusted H, referred to as H<sub>a</sub>, can be obtained by postulating that the increase in hits dH per unit increase of forecasts dF is proportional to the difference between observations and hits, that is, dH/dF ∝ (O - H). This relation proceeds from the self-evident statement that the larger the area (# of grid boxes) covered by misses (O - H), the larger is the probability that an additional 'yes' forecast provides an additional hit. Nevertheless, this formulation may lead to the paradoxical result that H > F (Mesinger, 2008); thus, as an alternative postulation, the increase in hits dH per unit increase of false alarms dA (= dF - dH) can be set as proportional to the unhit area, that is, dH/dA ∝ (O - H).

More details on the verification measures used are available in Wilks (2006), Jolliffe and Stephenson (2011), WWRP/WGNE Joint Working Group on Forecast Verification Research (2013) and the aforementioned references.

## 4.2. Bootstrap-based confidence intervals

To evaluate whether differences in scores between the two competing model versions were statistically significant, the bootstrap-based hypothesis test introduced by Hamill (1999) was used. Commonly applied hypothesis tests usually require conditions concerning the probability distribution of differences in score that are rarely satisfied. The basic idea of the bootstrap-based techniques is to build an artificial dataset by repeating random sampling, with replacement, from the original sample data, without making any assumption on the underlying probability distribution (Dionnis and Efron, 1983). However, this technique, like other hypothesis test methods, requires that the sample data considered have negligible autocorrelation in space and time. In the present case, the contingency table elements used to compute the categorical scores were tallied over the whole set of grid points available each day in Italy for the purpose of preventing the spatial autocorrelation affecting the hypothesis test results (see Hamill, 1999). Furthermore, score time series were negligibly temporally autocorrelated when observations and forecasts were accumulated on a daily basis (not shown).

The bootstrap-based test was then applied as follows. Let S be one of the scores calculated, and the null hypothesis H<sub>0</sub> of the Hamill test is that S<sub>M1</sub> - S<sub>M2</sub> = 0, where the subscripts M1 and M2 indicate the two competing NWP models, namely QBOLAM-KF and BOLAM-KF. The alternative hypothesis H<sub>A</sub> is given by S<sub>M1</sub> - S<sub>M2</sub> ≠ 0. The bootstrap method is then performed by randomly generating an indexed family {I<sub>t</sub>}<sub>t=1,...,T</sub> over the entire set of considered T days (= 183 days in this study), where each element can be equal to either 1 or 2. For a given day t, the two daily contingency tables for M1 and M2 are swapped if the index I<sub>t</sub> is equal to 2, whereas the swapping is not performed if I<sub>t</sub> is equal to 1. Using these random indicators, the resampled scores for the two competing configurations, S<sup>*</sup><sub>M1</sub> and S<sup>*</sup><sub>M2</sub>, are then calculated using the sum of the elements of the shuffled T daily contingency
tables. By repeating the procedure $N_B (=10^4)$ times, $N_B$ different resampled scores are built for each model version. The null hypothesis $H_0$ is then evaluated by means of a two-tailed test using the percentile method (Wilks, 2006) and assuming 0.05 as the $\alpha$ significance level. Therefore, the 95% confidence interval is determined by verifying the position of the score differences $S_{M1} - S_{M2}$ in the resampled distribution of $S_{M1}^* - S_{M2}^*$, that is, by defining the largest and the smallest $N_B \alpha/2$ of the $N_B$ bootstrap sample differences. Score differences outside the so-obtained confidence intervals are considered statistically significant for the chosen $\alpha$; in other words, the hypothesis $H_0$ is rejected.

In the present study, the confidence intervals are graphically referenced to QBOLAM-KF, because this version was selected as the benchmark reference. Nevertheless, results do not change when referencing the confidence intervals to the competitor BOLAM-KF forecasting system. The Hamill test is a transpose symmetric technique that is invariant under interchanging the roles of the two competing forecasting systems (Accadia et al., 2003a).

### 4.3. Fourier power spectrum analysis

As described in Section 4.1, differences in bias should be taken into account when comparing categorical scores to prevent any misinterpretation of the performance of the competing models. In addition to this, another complication derives from the differences in the spatial scales resolved by the NWP models under investigation. In the literature, several studies (e.g. Mass et al., 2002; Weygandt et al., 2004; Gilleland, 2013) have pointed out that detailed forecast fields that contain differing degrees of small-scale detail could score worse than smoother forecast fields because of the so-called double penalty effect, due to sensitivity to small displacement errors. In addition, NWP models tend to under-represent the smallest scales, resulting in a real resolution coarser than the nominal grid mesh size (see, e.g. Grasso, 2000; Harris et al., 2001; Chèruy et al., 2004; Frehlich and Sharman, 2008).

A similar result was found in another study when comparing the precipitation forecasts produced by three different NWP models (Lanciani et al., 2008). Even if the three models had the same nominal grid spacing, it was demonstrated that they differed in terms of both the real spatial resolution and the amount of small-scale details, and such differences had a remarkable impact on verification measures. A reliable intercomparison was instead achieved by up-scaling (see also Cherubini et al., 2002) the forecast fields on a coarser verification grid; in this way, differences among the up-scaled forecasts' spectra were smaller and negligible in terms of impact on scores.

In order to assess if the two competing BOLAM-based forecasting systems were defined on grids with the same resolution and if they had the same amount of small-scale details, the suggestion of Weygandt et al. (2004) to perform a Fourier power spectrum analysis was adopted.

The spectrum analysis can be performed as follows. Given a daily precipitation field $\Phi(x,y)$, the corresponding 2D energy spectrum $E(k_x, k_y)$ is then defined formally in terms of the $k_x$ and $k_y$ wavenumber components as the Fourier transform of the autocorrelation function $K(l_x, l_y)$ at lag $l_x$ in the x direction and at lag $l_y$ in the y direction (Harris et al., 2001; Press et al., 2007). However, according to the Wiener-Khinchin theorem, the Fourier power spectrum can be computed by multiplying

the 2D Fourier transform by its complex conjugate, that is:

$$E(k_x, k_y) = \frac{1}{2\pi} \left| \int e^{-i(k_x x + k_y y)} \Phi(x, y) \, dx \, dy \right|^2 \quad (1)$$

This calculation is preferred, because it suppresses some computational noise. In addition, a Hann window should be applied to filter forecast data and reduce aliasing.

The 2D spectrum can then be presented as an isotropic power spectrum $E(k)$, where $k = \sqrt{k_x^2 + k_y^2}$, through angular averaging (Equation 1). Scaling of the power spectrum occurs when $E(k) = k^{-\beta}$. In other words, the spectrum shows scale invariance if it is linear in $k$ on a log-log plot. The spectral slope $\beta$ is then an indicator of the forecast field's smoothness: the higher the $\beta$ is, the smoother, more organized the structure of the forecast field is, with the spectral energy mainly concentrated in the mean and large scales (Davis et al., 1996; Harris et al., 1996). According to the recent literature (Koh et al., 2012), it is not possible to assign a unique range of variability for the empirical $\beta$ values associated with the precipitation fields. The spectral scaling exponent appears to be highly dependent on several factors, such as the accumulation time step considered, the regional climate of the selected verification area, and the rain-generating mechanism.

In addition, it may become unfeasible to evaluate the single spectral scaling exponent when verifying a long set of forecast fields. Therefore, a mean scaling slope $\bar{\beta}$ ($\pm$standard deviation) was introduced to characterize the structure of the QPFs under investigation by averaging the daily $\beta$ values over the verification time period (April to September 2001). This is somewhat similar to that done by Koh et al. (2012) when validating the QPFs at different gauge locations for a single event.

### 4.4. The 'bubble' approach for the spatial distribution of the contingency table frequencies

The contingency table frequencies reported in Table 2 were also used to assess qualitatively whether the forecast performance of the two compared BOLAM model versions was dependent on a particular location.

For selected thresholds, this is achieved by comparing the distributions of hits, false alarms and misses related to the two forecast series when clustered at different verification grid points (Accadia et al., 2005; Mariani et al., 2005). Graphically, a map is produced by plotting in each grid point a 'bubble' whose size is proportional to the number of occurrences of the element of the contingency table under investigation. Considering that the correct no-rain forecasts are usually distributed almost uniformly over the verification grid, its spatial distribution is not considered in the 'bubble' approach, because it does not provide any particular information on the performance of a forecasting system.

## 5. Results of the intercomparison

As the first phase of the verification study QPFs modelled by BOLAM-KF and QBOLAM-KF were evaluated against the corresponding observed rainfall analyses over the native 0.1° grid domain. The verification scores reported in Table 3 were computed from the sum of the 183 daily contingency tables to reduce their sensitivity to small changes in the contingency table population. Differences in scores were tested for statistical
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 4. The 24 h accumulation scores, reported as a function of the selected thresholds, for a set of 183 days over Italy for BOLAM-KF ('competitor' – solid line) and QBOLAM-KF ('reference' – dashed line) on the native 0.1° model grid: (a) BIAS, (b) ETS, (c) HK, (d) POD, and (e) FAR. The confidence intervals are referenced to QBOLAM-KF scores and they indicate the 2.5th and the 97.5th percentiles of the bootstrap-based distributions of the differences in scores.</figcaption></figure>

significance by applying the Hamill bootstrap-based hypothesis test with a 5% significance level. Figures 4–7 illustrate the score comparisons. In addition, the frequency bias values obtained for the two model versions with and without the use of the BIAS adjustment method are reported in Table 4.

The initial focus is on the traditional scores calculated with neither recalibration nor BIAS adjustment (Figure 4). In Figure 4(a), the updated BOLAM version has a BIAS close to 1 for almost all the considered thresholds. As depicted in Table 4, its BIAS values range from a minimum of 0.879 (for 40.0 mm 24 h⁻¹) to a maximum of 1.355 (for 0.5 mm 24 h⁻¹). On the contrary, the QBOLAM version shows a BIAS score above 1 for all the thresholds, reaching a maximum value equal to about 2.7 at the highest thresholds (Table 4), which means that this model version tends to overpredict the occurrence of rainfall events at any rate. Moreover, differences in BIAS between these two versions turn out to be always statistically significant, except for 0.5 mm 24 h⁻¹. As mentioned in Section 2, the former QBOLAM-Kuo version proved to be even worse

than QBOLAM-KF in the same summertime, especially when considering the higher rainfall thresholds, reaching a maximum BIAS values equal to 5.4 (see BIAS results in Mariani, 2008).

In terms of ETS and FAR (Figures 4(b) and (e)) BOLAM-KF scores better than QBOLAM-KF and the bootstrap-based confidence intervals show that the existing differences between the two model versions are statistically significant for all the thresholds. This is not the case for HK and POD. For HK, BOLAM-KF appears to perform statistically significantly better than the QBOLAM-KF only for the lower thresholds (Figure 4(c)). For the higher thresholds, the HK values for the two model versions are close to each other, seeming to indicate that there is no difference in the forecast performance when predicting low probability rain events. For POD, the differences in the model versions' performance appear not to be statistically significant at all, apart from the lowest rainfall threshold (Figure 4(d)).

The interpretation of the verification results is clearer when considering the EDI and SEDI scores, which were calculated
<figure><img src="image_12.png" /><figcaption>Figure 5. As in Figure 4, but using recalibrated QPFs for BOLAM-KF (solid line) and QBOLAM-KF (dashed line): (a) EDI and (b) SEDI.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 6. As in Figure 4, but after applying the BIAS adjustment method to give frequency bias values of QBOLAM-KF (dashed line) similar to those of BOLAM-KF (solid line): (a) ETS, (b) HK, (c) POD, and (d) FAR.</figcaption></figure>

over the recalibrated forecasts (→BIAS = 1). Figure 5 shows that BOLAM-KF performs better than QBOLAM-KF for all the thresholds and, in particular, for the thresholds associated with the low-base rate rainfall events: the values for BOLAM-KF are at least 1.3 times greater than those computed for QBOLAM-KF. The differences in EDI and SEDI between the two competing forecasting systems appear not to be statistically significant only for the highest threshold where, however, they are near to the limits of the bootstrap-based confidence intervals.

On the basis of these results, the effect of the differences in bias on the traditional scores was then studied, in particular, on HK and POD that showed a different picture of the performance capability of the two model versions. Two sets of scores were considered: (1) the ETS, HK, POD and FAR (see Figure 6 and Table 4) were recalculated after relabelling the forecast contours of the more biased model version, QBOLAM-KF, so as to obtain BIAS values similar to those of the BOLAM-KF

model version (this was done via the Hamill BIAS adjustment methods); and (2) the bias-adjusted $ETSA_{dH/dF}$ and $ETSA_{dH/dA}$ (see Figure 7) that by definition remove the bias dependency. In both cases, the scores prove undoubtedly that the QBOLAM-KF version scores are worse than the BOLAM-KF version and the differences in scores are statistically significant for all the rainfall thresholds, except for HK and POD. For these measures, the differences at the higher thresholds are not statistically significant even though they are very near to the limits of the confidence intervals (Figures 6(b) and (c)).

It was investigated whether the differences in scores were in some way affected by the differences in the amount of scales represented in each of the two QPF series. Therefore, the mean scaling spectral slopes of the precipitation forecast fields for QBOLAM-KF (hereinafter $\bar{\beta}_{QBOLAM}$) and BOLAM-KF (hereinafter $\bar{\beta}_{BOLAM}$) were calculated by averaging the corresponding daily $\beta$ values over the verification time period, which were achieved by scaling the isotropic discrete power
<figure><img src="image_14.png" /><figcaption>Figure 7. The 24 h accumulation bias-adjusted ETS scores for a set of 183 days over Italy for BOLAM-KF (solid line) and QBOLAM-KF (dashed line) on the native 0.1° model grid: (a) ETSAdH/DF and (b) ETSAdH/DA. The confidence intervals are referenced to QBOLAM-KF scores and they indicate the 2.5th and the 97.5th percentiles of the bootstrap-based distributions of the differences in scores.</figcaption></figure>

<table><caption>Table 4. Comparison of the frequency bias values between BOLAM-KF and QBOLAM-KF computed over both the native 0.1° model grid and the 0.5° verification grid, without and with using the Hamill BIAS adjustment method.</caption><thead><tr><th>Grid size</th><th>Original threshold (mm 24 h<sup>-1</sup>)</th><th>BIAS<sup>a</sup> of BOLAM-KF</th><th>BIAS<sup>a</sup> of QBOLAM-KF</th><th>BIAS difference</th><th>Adjusted threshold (mm 24 h<sup>-1</sup>)</th><th>Adjusted BIAS<sup>b</sup> of QBOLAM-KF</th><th>Adjusted BIAS difference</th></tr></thead><tbody><tr><td rowspan="6">0.1°</td><td>0.5</td><td>1.355</td><td>1.365</td><td>-0.010</td><td>0.5<sup>c</sup></td><td>1.365</td><td>-0.010</td></tr><tr><td>5.0</td><td>1.193</td><td>1.780</td><td>-0.587</td><td>8.9</td><td>1.243</td><td>-0.050</td></tr><tr><td>10.0</td><td>1.158</td><td>2.034</td><td>-0.876</td><td>16.5</td><td>1.202</td><td>-0.044</td></tr><tr><td>20.0</td><td>1.060</td><td>2.420</td><td>-1.360</td><td>30.5</td><td>1.097</td><td>-0.037</td></tr><tr><td>30.0</td><td>1.046</td><td>2.755</td><td>-1.709</td><td>44.0</td><td>1.079</td><td>-0.033</td></tr><tr><td>40.0</td><td>0.879</td><td>2.679</td><td>-1.800</td><td>58.0</td><td>0.912</td><td>-0.033</td></tr><tr><td rowspan="6">0.5°</td><td>0.5</td><td>1.352</td><td>1.720</td><td>-0.368</td><td>1.3</td><td>1.379</td><td>-0.027</td></tr><tr><td>5.0</td><td>1.198</td><td>1.819</td><td>-0.621</td><td>8.1</td><td>1.198</td><td>-0.025</td></tr><tr><td>10.0</td><td>1.154</td><td>1.863</td><td>-0.709</td><td>14.5</td><td>1.154</td><td>0.011</td></tr><tr><td>20.0</td><td>1.058</td><td>1.906</td><td>-0.848</td><td>26.0</td><td>1.058</td><td>-0.010</td></tr><tr><td>30.0</td><td>0.989</td><td>1.906</td><td>-0.917</td><td>37.5</td><td>0.989</td><td>0.011</td></tr><tr><td>40.0</td><td>0.663</td><td>1.791</td><td>-1.128</td><td>51.0</td><td>0.663</td><td>-0.046</td></tr></tbody></table>

spectrum $E(k)$ in Equation 1 using the linear least-squares fitting technique. The least-squares fit of the first part of the power spectrum was applied considering different numbers of fitted points, namely the first 5, 7 and 12 points to exclude any dependence of this parameter on the results. This implied that a spatial window was considered for the spectral analysis spanning from 707.1 km to 78.6, 54.4 and 30.7 km, respectively. For each number of fitted points, the goodness of the fit was measured by computing the mean of the root mean square errors (RMSEs) associated with the least-squares fits.

Figure 8 reports the comparison of $\bar{\beta}_{QBOLAM}$ versus $\bar{\beta}_{BOLAM}$. For completeness the mean scaling slopes of the QPFs modelled by QBOLAM-Kuo are provided as well. Apart from the variability of the scaling slopes (see standard deviation intervals in Figure 8) that was caused by the different weather systems that occurred during the 6 months, a strong average difference is noticeable in the structure of QPFs among the compared forecasting systems. The mean scaling slope ranges from 1.97 to 2.61 for BOLAM-KF and from 1.55 to 2.06 for both QBOLAM versions, depending on the number of fitted points. Thus, on an average, the precipitation fields generated with the latter operational model version are smoother and more organized over the native grid domain than those modelled by the former QBOLAM versions. A similar outcome is obtained when scaling the

second part of the power spectrum (not shown). In other words, the amount of scale details present in the BOLAM-KF QPFs differ from those present in the QPFs generated with the two QBOLAM versions when considering the native 0.1° grid as verification grid domain. As previously pointed out, intercomparing forecasts with different power spectra could lead to a misrepresentation of the actual performance of the forecasting systems investigated. In order to remove doubt that BOLAM-KF scored better because of its general tendency to generate smoother, more organized precipitation fields over the native 0.1° grid (i.e. $\bar{\beta}_{QBOLAM} > \bar{\beta}_{QBOLAM}$), it was decided to repeat the intercomparison using a coarser verification grid selected in such a way so as to guarantee the similarity of the mean scaling slopes (i.e. $\bar{\beta}_{QBOLAM} \approx \bar{\beta}_{BOLAM}$). This was achieved by remapping (Accadia et al., 2003b) the 24 h precipitation forecast time series modelled by QBOLAM-KF and BOLAM-KF onto a coarser 0.5° grid, because for both versions these remapped QPF series had a mean scaling slope of about 3.7 (not shown).

The categorical scores were hence calculated over the 0.5° verification grid (see Figures 9–12 and Table 4). The outcomes, especially those obtained after removing or addressing the BIAS effect (Figures 10–12), confirm what has been described previously for the intercomparison over the native 0.1° grid: BOLAM-KF scores better than QBOLAM-KF in predicting

<footnote><sup>a</sup>The BIAS score for the two competing models was originally computed by using the same rainfall threshold for both observations and forecasts.</footnote>

<footnote><sup>b</sup>The adjusted BIAS score for QBOLAM-KF was obtained by introducing an adjusted threshold for forecasts, while keeping unchanged the original threshold for observations.</footnote>

<footnote><sup>c</sup>Forecast threshold not adjusted, because the absolute difference between the BIAS values of the two competing models was already less than the error term $\epsilon_b$ equal to 0.05.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 8. Comparison of the mean scaling slope $\bar{\beta}$ for the 24 h QPFs modelled by QBOLAM-KF and BOLAM-KF computed by averaging the corresponding daily $\beta$ values over April to September 2001, which were obtained by fitting the first part of the Fourier power spectra $E(k)$ with respect to the first 5, 7, and 12 fitted points. The standard deviation bars of the $\bar{\beta}$ values and the mean values of the daily RMSEs related to the least-squares fits are reported. For completeness, the mean scaling slopes for the QPFs generated with QBOLAM-Kuo are presented as well.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 9. As in Figure 4, but for the 24 h accumulation scores calculated for a set of 183 days over Italy for BOLAM-KF (solid line) and QBOLAM-KF (dashed line) remapped on the 0.5° verification grid: (a) BIAS, (b) ETS, (c) HK, (d) POD, and (e) FAR.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 10. As in Figure 9, but using recalibrated 0.5° QPFs for BOLAM-KF (solid line) and QBOLAM-KF (dashed line): (a) EDI and (b) SEDI.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 11. As in Figure 9, but after applying the BIAS adjustment method to give the frequency bias values of QBOLAM-KF (dashed line) similar to those of BOLAM-KF (solid line) over the 0.5° verification grid: (a) ETS, (b) HK, (c) POD, and (d) FAR.</figcaption></figure>

both high-base rate and low-base rate precipitation events. The former operational version is penalized by its tendency to produce a high number of false alarms (see Figures 9(a) and (e) and the frequency bias values in Table 4).

Except for the bias-adjusted ETS scores that show that the differences between the updated model version and the older model version are statistically significant for each one of the considered rainfall thresholds (Figure 12), the other differences in scores are statistically significant only at the low and medium thresholds (Figures 10 and 11). The differences at the higher thresholds are instead very near to the boundaries of the confidence intervals, so that they cannot be considered statistically significant. Nevertheless, the wide confidence intervals obtained at the higher thresholds, which cause the acceptance of the null hypothesis *H*<sub>0</sub>, might be a consequence of an inadequate (small) sample size of the rare precipitation events that affect, in particular, the 0.5° intercomparison. This issue could be investigated further in the future by considering a longer dataset of observed and forecast fields so as to guarantee an adequate sample size also for the rare events.

Additional information on the temporal variation in the forecast performances was obtained by computing the categorical scores over the single-month datasets, using the 5.0 mm 24 h<sup>-1</sup> threshold to ensure robust samples for score calculation (Figure 13). The monthly frequency biases for the 0.5° remapped QBOLAM-KF turn out to be always statistically greater than those of the 0.5° remapped BOLAM-KF, with two large peaks (BIAS ≅ 3) in July and August 2001 (Figure 13(a)). These within-season differences confirm the QBOLAM-KF tendency to overpredict the occurrence of rain events in the summer months. This is also valid when evaluating the two model versions over the native 0.1° grid (not shown). Even when the effect on score comparison of the monthly differences in BIAS is removed *via* either the use of the bias-adjusted threat score (Figure 13(b)) or the application of the Hamill BIAS adjustment procedure (Figures 13(c) and (d)), the QBOLAM-KF version continues to perform worse than the BOLAM-KF version. Differences in score, which are on an average statistically significant or close to boundaries of the bootstrap-based intervals, show a decrease in the forecast performance of
<figure><img src="image_19.png" /><figcaption>Figure 12. As in Figure 7, but for the bias-adjusted ETS scores calculated for a set of 183 days over Italy for BOLAM-KF (solid line) and QBOLAM-KF (dashed line) on the 0.5° verification grid: (a) $ETSA_{dH/dF}$ and (b) $ETSA_{dH/dA}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 13. The monthly 24 h accumulation scores over Italy for BOLAM-KF (solid line) and QBOLAM-KF (dashed line) remapped on the 0.5° grid, computed using the 5.0 mm 24 h⁻¹ threshold: (a) BIA, (b) $ETSA_{dH/dA}$, (c) HK, and (d) FAR. The confidence intervals are referenced to QBOLAM-KF scores and they indicate the 2.5th and the 97.5th percentiles of the bootstrap-based distributions of the differences in scores. The QBOLAM-KF scores in panels (c) and (d) are calculated after applying the Hamill BIAS adjustment against BOLAM-KF.</figcaption></figure>

QBOLAM-KF between June and August, which is associated with a general increase in false alarms (Figures 13(b)–(d)). It is worth noting that BOLAM-KF-localized increase (decrease) in FAR in July (September) has, as a consequence, a decrease (increase) in the $ETSA_{dH/dA}$ value compared with the general behaviour of the skill score (Figures 13(d) vs. (b)); this is owing to a relative increase in the number of false alarms (hits). By contrast, the observed increase in the June value of the BOLAM-KF HK is related to the relative decrease of misses as depicted by the POD month-to-month variation (not shown).

Finally, the ‘bubble’ approach was applied to assess how the two model versions spatially distribute the 0.1° and 0.5° grid-based contingency table frequencies related to the 5.0, 10.0, and 20.0 mm 24 h⁻¹ rainfall thresholds. Except from the differences in the results associated with the different grids and thresholds, the ‘bubble’ approach clearly and visually identifies the same areas of improvement of the new BOLAM version. A comparison between the distributions of hits, false alarms, and misses for the 0.5° remapped QBOLAM-KF QPFs and the corresponding distributions for the 0.5° remapped BOLAM-KF

QPFs for a rainfall threshold of 10.0 mm 24 h⁻¹ is analyzed (see Figure 14). For BOLAM-KF, a large decrease in false alarms is observed especially over the Italian mountain ranges (Figures 14(c) vs. (d)) and a slight general decrease in misses (Figures 14(e) vs. (f)). No particular change is identified in the spatial distribution of hits (Figures 14(a) vs. (b)), suggesting that during the analysed period most of the correct no-rain predictions were forecasted erroneously as events by the previously operational QBOLAM version (not shown).

## 6. Conclusions and final remarks

This paper discusses the quantitative precipitation forecast (QPF) performance of the new version of the meteorological model BOlogna Limited Area Model (BOLAM) that deploys more efficient and advanced schemes for the parameterization of cumulus convection, radiation, soil, and exchange of turbulent fluxes, and an explicit microphysical scheme with five hydrometeors. This up-to-date model version is currently operational within the Institute for Environmental Protection and
<figure><img src="image_21.png" /><figcaption>Figure 14. Hits, false alarms and misses for QBOLAM-KF (left column) and BOLAM-KF (right column) calculated over the 0.5° verification grid for the period from 1 April to 30 September 2001 by using the 10.0 mm 24 h⁻¹ rainfall threshold.</figcaption></figure>

Research (ISPRA) Hydro-Meteo-Marine Forecasting System (Sistema Idro-Meteo-Mare – SIMM) forecasting system.

A thorough forecast verification was accomplished by using a multimethod approach based on: (1) a power spectrum analysis to assess the version differences in terms of spatial scale details; (2) an intercomparison of several categorical scores to measure the model performance when dealing with both frequent and rare precipitation events; (3) bootstrap-based confidence intervals to evaluate the statistical significance of the

differences in the score performance; and (4) a 'bubble' chart analysis to display the spatial distribution of the contingency table elements. The verification outcomes for the new BOLAM version have been easily connected to those related to the previous model versions by considering a past verification time period: April to September 2001. This period was the most critical for the performance of the two Quadrics BOlogna Limited Area Model (QBOLAM) versions owing to an excessive occurrence of precipitation false alarms. As a consequence,
these previous versions were characterized by frequency bias values much greater than 1. On the contrary, in this study, the currently operational BOLAM version has exhibited a frequency bias almost equal to 1 over the same time period. Such a difference can affect the score intercomparison providing ambiguous results of the actual performance capabilities of the compared model versions, unless the bias dependency would be somehow considered as it was done in this statistical analysis.

The comparison of the QPF performances of BOLAM-KF versus QBOLAM-KF was tested over two different grids: the native 0.1° model grid and a coarser 0.5° grid, which was selected so as to prevent, similarly, the verification results from being affected by the differences in the scales resolved by the two model versions. In both cases, the BOLAM-KF version shows a statistically significantly better performance in predicting both high-base rate and low-base rate precipitation events. An analysis of the month-to-month variation of the categorical scores confirms this picture highlighting, in addition, a large deficiency for the QBOLAM-KF performance in July and August.

On the basis of these positive outcomes, some of the authors of the present study have decided to examine the potential of increasing the BOLAM forecast quality by improving the model configuration, in particular, in terms of horizontal grid spacing, domain extension, data quality and resolution of initial and boundary conditions (IC/BCs), and nesting design. A first series of experiments has been conducted recently to determine whether these modifications lead to improvements (Casaioli et al., 2013). The results have highlighted the effectiveness of decreasing the BOLAM grid length to 0.07° and simultaneously increasing the domain extension to improve the QPF quality, provided that higher resolution IC/BCs are applied directly to the 0.07° BOLAM without using any coarser resolution parent model run. Future experiments will examine separately the effect of increasing the domain extension and the effect of increasing model resolution. In addition, it will be assessed how the performance results change as a function of lead time.

Furthermore, future work will investigate by means of feature-based methods, how the new BOLAM version accounts for the spatial structure of weather forecast fields. In particular, the idea is to consider different verification time periods related to international monitoring and forecasting initiatives on high impact weather events, like those endorsed by the World Meteorological Organization. Such initiatives can indeed ensure the availability of measurements over larger domains and the accessibility to both traditional (e.g. rain gauge) and new sources of data (e.g. radar). This can also be helpful in overcoming shortcomings, such as the well-documented 'domain jumping' (Tartaglione et al., 2005; Grams et al., 2006; Mariani et al., 2008; Ebert and Gallus, 2009), which mainly occurs when spatial verification methods are applied to small or narrow domains, as in the present case of the Italian peninsula.

## Acknowledgements

This work was performed under the framework of the SIMM forecast verification programme and it includes some of the research undertaken by Dr Stefano Mariani during his PhD studies at the University of Ferrara, Italy. We would like to thank the Regional MET Services of Calabria, Emilia-Romagna, Ligure, Marche, Piedmont, Sardinia, Sicily, and Valle d'Aosta for providing additional rain gauge data. A special thanks goes to Dr Piero Malguzzi (ISAC-CNR, Bologna, Italy) and

Prof. Mariarosaria Padula (University of Ferrara, Italy), who supported the research activity with fruitful suggestions, and to Mr Stefano Tagliaventi (CINECA, Rome, Italy), who collaborated on the updating of the SIMM forecasting system that was conducted by ISPRA together with ISAC-CNR. Finally, we gratefully acknowledge helpful comments and suggestions that we received from two anonymous reviewers of this journal.

## References

Accadia C, Casaioli M, Mariani S, Lavagnini A, Speranza A, Di Venere A, Inghilesi R, Ferretti R, Paolucci T, Cesari D, Patruno P, Boni G, Bovo S, Cremonini R. 2003a. Application of a statistical methodology for limited area model intercomparison using a bootstrap technique. *Nuovo Cimento* **26C**: 61–77.

Accadia C, Mariani S, Casaioli M, Lavagnini A, Speranza A. 2003b. Sensitivity of precipitation forecast skill scores to bilinear interpolation and a simple nearest-neighbor average method on high-resolution verification grids. *Weather Forecast*. **18**: 918–932.

Accadia C, Mariani S, Casaioli M, Lavagnini A, Speranza A. 2005. Verification of precipitation forecasts from two limited-area models over Italy and comparison with ECMWF forecasts using a resampling technique. *Weather Forecast*. **20**: 276–300.

Ahrens B, Beck A. 2008. On upscaling of rain-gauge data for evaluating numerical weather forecasts. *Meteorol. Atmos. Phys*. **99**: 155–167.

Barnes SL. 1964. A technique for maximizing details in numerical weather map analysis. *J. Appl. Meteorol.* **3**: 396–409.

Barnes SL. 1973. Mesoscale objective analysis using weighted time-series observations. NOAA Technical Memorandum ERL NSSL-62, Environmental Research Laboratories, National Severe Storms Laboratory: Norman, OK; 60.

Blumberg AF, Mellor GL. 1987. A description of a three-dimensional coastal ocean circulation model. In *Three-Dimensional Coastal Ocean Models*, Coastal and Estuarine Science Series, Heaps NS (ed), Vol. 4. American Geophysical Union: Washington, DC; 1–16.

Brooks HE, Doswell CA III. 1993. New technology and numerical weather prediction – A wasted opportunity? *Weather* **48**: 173–177.

Buzzi A, Fantini M, Malguzzi P, Nerozzi F. 1994. Validation of a limited area model in cases of Mediterranean cyclogenesis: surface fields and precipitation scores. *Meteorol. Atmos. Phys*. **53**: 53–67.

Casaioli M, Mariani S, Malguzzi P, Speranza A. 2013. Factors affecting the quality of QPF: a multi-method verification of multi-configuration BOLAM reforecasts against MAP D-PHASE observations. *Meteorol. Appl*. **20**: 150–163.

Casati B, Wilson LJ, Stephenson DB, Nurmi P, Ghelli A, Pocernich M, Damrath U, Ebert EE, Brown BG, Mason S. 2008. Forecast verification: current status and future directions. *Meteorol. Appl*. **15**: 3–18.

Cherubini T, Ghelli A, Lalaurette F. 2002. Verification of precipitation forecasts over the Alpine Region using a high-density observing network. *Weather Forecast*. **17**: 238–249.

Chêruy F, Speranza A, Sutera A, Tartaglione N. 2004. Surface winds in the Euro-Mediterranean area: the real resolution of numerical grids. *Ann. Geophys*. **22**: 4043–4048.

Davis A, Marshak A, Wiscombe W, Cahalan R. 1996. Scale invariance of liquid water distributions in marine stratocumulus. Part I: spectral properties and stationarity issues. *J. Atmos. Sci*. **53**: 1538–1558.

Davis C, Brown B, Bullock R. 2006a. Object-based verification of precipitation forecasts. Part I: methodology and application to mesoscale rain areas. *Mon. Weather Rev*. **134**: 1772–1784.

Davis C, Brown B, Bullock R. 2006b. Object-based verification of precipitation forecasts. Part II: application to convective rain systems. *Mon. Weather Rev*. **134**: 1785–1795.

Diaconis P, Efron B. 1983. Computer-intensive methods in statistics. *Sci. Am*. **248**: 116–130.

Ebert EE. 2001. Verification of precipitation areas. In *Proceedings of the 1st WWRP/WMO Workshop on the Verification of Quantitative Precipitation Forecasts*, 14–16 May 2001, Prague, Czech Republic. http://www.cawcr.gov.au/projects/verification/QPF_Prague2001/papers/beth_full_paper.doc (accessed 16 August 2013).

Ebert EE. 2008. Fuzzy verification of high-resolution gridded forecasts: a review and proposed framework. *Meteorol. Appl*. **15**: 51–64.
Ebert EE, Gallus WA. 2009. Toward better understanding of the Contiguous Rain Area (CRA) method for spatial forecast verification. *Weather Forecast*. **24**: 1401–1415.

Ebert EE, McBride JL. 2000. Verification of precipitation in weather systems: determination of systematic errors. *J. Hydrol*. **239**: 179–202.

Ebert EE, Damrath U, Wergen W, Baldwin ME. 2003. The WGNE assessment of short-term quantitative precipitation forecasts. *Bull. Am. Meteorol. Soc*. **84**: 481–492.

Ferro CAT, Stephenson DB. 2011. Extremal dependence indices: improved verification measures for deterministic forecasts of rare binary events. *Weather Forecast*. **26**: 699–713.

Frehlich R, Sharman R. 2008. The use of structure functions and spectra from numerical model output to determine effective model resolution. *Mon. Weather Rev*. **136**: 1537–1553.

Gallus WA Jr, Segal M. 2001. Impact of improved initialization of mesoscale features on convective system rainfall in 10-km Eta simulations. *Weather Forecast*. **16**: 680–696.

Ghelli A, Lalaurette F. 2000. Verifying precipitation forecasts using upscaled observations. *ECMWF Newsl*. **87**: 9–17.

Gilleland E. 2013. Testing competing precipitation forecasts accurately and efficiently: the spatial prediction comparison test. *Mon. Weather Rev*. **141**: 340–355.

Gilleland E, Ahijevych D, Brown BG, Casati B, Ebert EE. 2009. Intercomparison of spatial verification methods. *Weather Forecast*. **24**: 1416–1430.

Gilleland E, Ahijevych D, Brown BG, Ebert EE. 2010. Verifying forecasts spatially. *Bull. Am. Meteorol. Soc*. **91**: 1365–1373.

Göber M, Zsótér E, Richardson DS. 2008. Could a perfect model ever satisfy the forecaster? On grid box mean versus point verification. *Meteorol. Appl*. **15**: 359–365.

Grams JS, Gallus WA, Koch SE, Wharton LS, Loughe AF, Ebert EE. 2006. The use of a modified Ebert-McBride technique to evaluate mesoscale model QPF as a function of convective system morphology during IHOP 2002. *Weather Forecast*. **21**: 288–306.

Grasso LD. 2000. The difference between grid spacing and resolution and their application to numerical modeling. *Bull. Am. Meteorol. Soc*. **81**: 579–580.

Hagedorn R, Buizza R, Hamill TM, Leutbecher M, Palmer TN. 2012. Comparing TIGGE multimodel forecasts with reforecast-calibrated ECMWF ensemble forecasts. *Q. J. R. Meteorol. Soc*. **138**: 1814–1827.

Hamill TM. 1999. Hypothesis tests for evaluating numerical precipitation forecasts. *Weather Forecast*. **14**: 155–167.

Hamill TM, Whitaker JS, Wei X. 2004. Ensemble reforecasting: improving medium-range forecast skill using retrospective forecasts. *Mon. Weather Rev*. **132**: 1434–1447.

Hamill TM, Whitaker JS, Mullen SL. 2006. Reforecasts: an important dataset for improving weather predictions. *Bull. Am. Meteorol. Soc*. **87**: 33–46.

Hamill TM, Bates GT, Whitaker JS, Murray DR, Fiorino M, Galarneau TJ Jr, Zhu Y, Lapenta W. 2013. NOAA's second-generation global medium-range ensemble reforecast data set. *Bull. Am. Meteorol. Soc*. **94**: 1553–1565, DOI: 10.1175/BAMS-D-12-00014.1.

Harris D, Menabde M, Seed A, Austin GL. 1996. Multifractal characterization of rain fields with a strong orographic influence. *J. Geophys. Res*. **101**: 26405–26414.

Harris D, Foufoula-Georgiou E, Droegemeier KK, Levit JJ. 2001. Multiscale statistical properties of a high-resolution precipitation forecast. *J. Hydrometeorol*. **2**: 406–418.

Hoffman RN, Liu Z, Louis J-F, Grassotti C. 1995. Distortion representation of forecast errors. *Mon. Weather Rev*. **123**: 2758–2770.

Inghilesi R, Catini F, Bellotti G, Franco L, Orasi A, Corsini S. 2012. Implementation and validation of a coastal forecasting system for wind waves in the Mediterranean Sea. *Nat. Hazards Earth Syst. Sci*. **12**: 485–494.

Jolliffe IT, Stephenson DB. 2011. *Forecast Verification: A Practitioner's Guide in Atmospheric Science*, 2nd edn. John Wiley & Sons, Ltd.: Chichester; 292.

Kain JS, Fritsch JM. 1990. A one-dimensional entraining/detraining plume model and its application in convective parameterization. *J. Atmos. Sci*. **47**: 2784–2802.

Kain JS, Fritsch JM. 1993. Convective parameterization for mesoscale models: the Kain-Fritsch scheme. In *The Representation of Cumulus Convection in Numerical Models*, Raymond DJ, Emanuel KA (eds). American Meteorological Society: Boston, MA; 165–170.

Koch SE, desJardins M, Kocin PJ. 1983. An interactive Barnes objective map analysis scheme for use with satellite and conventional data. *J. Clim. Appl. Meteorol*. **22**: 1487–1503.

Koh T-Y, Bhatt BC, Cheung KKW, Teo CK, Lee YH, Roth M, Purwawirman. 2012. Using the spectral scaling exponent for validation of quantitative precipitation forecasts. *Meteorol. Atmos. Phys*. **115**: 35–45.

Kuo HL. 1974. Further studies of the parameterization of the influence of cumulus convection on large scale flow. *J. Atmos. Sci*. **31**: 1232–1240.

Lanciani A, Salvati M. 2008. Spatial interpolation of surface weather observations in Alpine meteorological service. FORALPS Technical Report 2, Università degli Studi di Trento, Dipartimento di Ingegneria Civile e Ambientale: Trento; 40. ISBN 978-88-8443-25-4. http://www.ing.unitt.it/dica/tools/download/Quaderni/Foralps_TR_2.pdf or http://opac.isprambiente.it/sebina/repository/catalogazione/immagini/pdf/Foralps_TR_2.pdf (accessed 25 January 2013).

Lanciani A, Mariani S, Casaioli M, Accadia C, Tartaglione N. 2008. A multiscale approach for precipitation verification applied to the FORALPS case studies. *Adv. Geosci*. **16**: 3–9.

Malguzzi P, Tartaglione N. 1999. An economical second-order advection scheme for numerical weather prediction. *Q. J. R. Meteorol. Soc*. **125**: 2291–2303.

Mariani S. 2008. Methodologies in research and operational contexts for numerical precipitation forecast verification, Doctoral thesis, University of Ferrara, Ferrara, Italy.

Mariani S, Casaioli M, Accadia C, Llasat MC, Pasi F, Davolio S, Elementi M, Ficca G, Romero R. 2005. A limited area model intercomparison on the "Montserrat-2000" flash-flood event using statistical and deterministic methods. *Nat. Hazards Earth Syst. Sci*. **5**: 565–581.

Mariani S, Accadia C, Tartaglione N, Casaioli M, Gabella M, Michaelides SC, Speranza A. 2008. Multi-sensor comparison and numerical modeling of atmospheric water fields: a VOLTAIRE case study over Cyprus. *Weather Forecast*. **23**: 674–701.

Mariani S, Casaioli M, Lanciani A, Tartaglione N, Accadia C. 2009. A multi-method intercomparison approach for precipitation fields modelled by LAMs in the Alpine area: two FORALPS case studies. *Meteorol. Atmos. Phys*. **103**: 79–92.

Mason I. 1989. Dependence of the critical success index on sample climate and threshold probability. *Aust. Meteorol. Mag*. **37**: 75–81.

Mass CF, Owens D, Westrick K, Colle AB. 2002. Does increasing horizontal resolution produce more skilful forecasts? *Bull. Am. Meteorol. Soc*. **83**: 407–430.

Mesinger F. 2008. Bias adjusted precipitation threat scores. *Adv. Geosci*. **16**: 137–142.

Mesinger F, Brill K. 2004. Bias normalized precipitation scores. In *Proceedings of 17th Conference on Probability and Statistics in the Atmospheric Sciences*, Seattle, WA, 11–15 January 2004. American Meteorological Society: Boston, MA.

Moise AF, Delage FP. 2011. New climate model metrics based on object-orientated pattern matching of rainfall. *J. Geophys. Rev*. **116**: D12108.

Murphy AH. 1993. What is a good forecast? An essay on the nature of goodness in weather forecasting. *Weather Forecast*. **8**: 281–293.

Nurmi P. 2003. Recommendations on the verification of local weather forecasts. ECMWF Technical Memorandum No. 430, ECMWF: Reading; 19.

Press WH, Teukolsky SA, Vetterling WT, Flannery BP. 2007. *Numerical Recipes: The Art of Scientific Computing*, 3rd edn. Cambridge University Press: Cambridge; 1256.

Raymond DJ, Emanuel KA. 1993. The Kuo cumulus parameterization. In *The Representation of Cumulus Convection in Numerical Models*, Raymond DJ, Emanuel KA (eds). American Meteorological Society: Boston, MA; 145–147.

Rossa A, Nurmi P, Ebert EE. 2008. Overview of methods for the verification of quantitative precipitation forecasts. In *Precipitation: Advances in Measurement, Estimation and Prediction*, Michaelides SC (ed). Springer-Verlag: New York, NY; 419–452.

Shuman FG. 1980. A modified threat score and a measure of placement error. NOAA/NWS Office Note 210; 13. http://www.ncep.noaa.gov/officenotes/NOAA-NPM-NCEPON-0003/014081E3.pdf (accessed 11 March 2013).

Skelly WC, Henderson-Sellers A. 1996. Grid box or grid point: what type of data do GCMs deliver to climate impacts researchers? *Int. J. Climatol*. **16**: 1079–1086.

Speranza A, Accadia C, Casaioli M, Mariani S, Monacelli G, Inghilesi R, Tartaglione N, Ruti PM, Carillo A, Bargagli A, Pisacane G, Valentinotti F, Lavagnini A. 2004. POSEIDON: an integrated system for analysis and forecast of hydrological, meteorological and surface marine fields in the Mediterranean area. *Nuovo Cimento 27C*: 329–345.
Speranza A, Accadia C, Mariani S, Casaioli M, Tartaglione N, Monacelli G, Ruti PM, Lavagnini A. 2007. SIMM: an integrated forecasting system for the Mediterranean area. *Meteorol. Appl.* **14**: 337–350.

Stephenson DB. 2000. Use of the “odds ratio” for diagnosing forecast skill. *Weather Forecast*. **15**: 221–232.

Tartaglione N, Mariani S, Accadia C, Speranza A, Casaioli M. 2005. Comparison of rain gauge observations with modeled precipitation over Cyprus using contiguous rain area analysis. *Atmos. Chem. Phys.* **5**: 2147–2154.

Tartaglione N, Mariani S, Casaioli M, Accadia C, Federicoe S, Michaelides SC. 2008. Searching for systematic location errors of quantitative precipitation forecasts over the Calabria region. *Meteorol. Appl.* **15**: 85–95.

The Wamdi Group. 1988. The WAM model-a third generation ocean waves prediction model. *J. Phys. Oceanogr.* **18**: 1776–1810.

Todini E, Ciarapica L. 2001. The TOPKAPI model. In *Mathematical Models of Large Watershed Hydrology*, Singh VP, Frevert DK, Meyer SP (eds). Water Resources Publications: Littleton, CO; 471–506.

Tustison B, Harris D, Foufoula-Georgiou E. 2001. Scale issues in verification of precipitation forecasts. *J. Geophys. Res.* **106**: 11775–11784.

Umgiesser G, Bergamasco A. 1993. A staggered grid finite element

model for the Venice Lagoon. In *Finite Elements in Fluids*, Morgan K, Ofiate E, Periaux J, Zienkiewicz OC (eds). Pineridge Press: Swansea; 659–668.

Umgiesser G, Melaku Canu D, Cucco A, Solidoro C. 2004. A finite element model for the Venice Lagoon. Development, set up, calibration and validation. *J. Mar. Syst.* **51**: 123–145.

Weygandt SS, Loughe AF, Benjamin SG, Mahoney JL. 2004. Scale sensitivities in model precipitation skill scores during IHOP. In *Proceedings of 22nd Conference on Severe Local Storms*, Hyannis, MA, 4–8 October 2004. American Meteorological Society: Boston, MA; 16A.8.

Wilks DS. 2006. *Statistical Methods in the Atmospheric Sciences*, International Geophysics Series, Vol. 91. Academic Press: San Diego, CA; 649.

WWRP/WGNE Joint Working Group on Forecast Verification Research. 2013. Forecast verification - issue, methods and FAQ. http://www.cawcr.gov.au/projects/verification/ (accessed 17 April 2013).

WWRP/WGNE Joint Working Group on Verification. 2004. Recommendations for the verification and intercomparison of QPFs from operational NWP models. http://www.cawcr.gov.au/projects/verification/WGNE/QPF_verif_recomm.pdf (accessed 20 August 2012).