# Calibrating the CAMS European multi-model air quality forecasts for regional
air pollution monitoring

Gabriele Casciaro<EMAIL_ADDRESS>Mattia Cavaiola
<EMAIL_ADDRESS>Andrea Mazzino<EMAIL_ADDRESS>
###### Abstract

The CAMS air quality multi-model forecasts have been assessed and calibrated
for PM10, PM2.5, O3, NO2, and CO against observations collected by the
Regional Monitoring Network of the Liguria region (northwestern Italy) in the
years 2019 and 2020. The calibration strategy used in the present work has its
roots in the well-established Ensemble Model Output Statistics (EMOS) through
which a raw ensemble forecast can be accurately transformed into a predictive
probability density function, with a simultaneous correction of biases and
dispersion errors. The strategy also provides a calibrated forecast of model
uncertainties. As a result of our analysis, the key role of pollutant real-
time observations to be ingested in the calibration strategy clearly emerge
especially in the shorter look-ahead forecast hours. Our dynamic calibration
strategy turns out to be superior with respect to its analogous where real-
time data are not taken into account. The best calibration strategy we have
identified makes the CAMS multi-model forecast system more reliable than other
raw air quality models running at higher spatial resolution which exploit more
detailed information from inventory emission. We expect positive impacts of
our research for identifying and set up reliable and economic air pollution
early warning systems.

###### keywords:

Ensemble forecasts , Air quality monitoring , Ensemble model output statistics
(EMOS) , Static vs. dynamic calibrations , Air pollution early warning systems

††journal: Atmospheric Environment

[DICCA]organization=DICCA, Department of Civil, Chemical and Environmental
Engineering. University of Genoa, addressline=Via Montallegro 1, city=Genoa,
postcode=16145, state=Genoa, country=Italy [INFN] organization=INFN, Nazional
Institute for Nuclear Physics, Genoa Section, addressline=Via Dodecaneso 33,
city=Genoa, postcode=16146, state=Genoa, country=Italy

A dynamic calibration strategy brings out the intrinsic strength of CAMS
multi-model forecasts

The CAMS calibrated forecasts neatly outperforms persistence

The CAMS calibrated forecasts has larger skills than raw higher-resolution
strategies

The CAMS calibrated system also provides calibrated forecasts of model
uncertainties

Our outcomes impact the management of air pollution early warning systems

##  1 Introduction

Among all environmental problems, air pollution is surely one of the most
serious causing, each year, at least 7 million premature deaths. There are
statistical evidences of positive correlations between deaths caused by air
pollution and the ground concentration of fine particulate matter, PM, (Yin et
al., 2017). Studies have also shown that PM, ozone, nitrogen dioxide, and
sulphur dioxide degrade ambient air quality and cause severe health problems
in humans (Yin et al., 2017). In way of example, increased ground-level ozone
is responsible for metabolic effects on humans leading to glucose intolerance
and hyperlipidemia (Shore, 2019).
Following the World Bank (2016), these dramatic consequences that affect
people’s lives also have huge implications for the world economies. The damage
to the local economies are estimated (World Bank, 2016) in US $ 5 trillion in
welfare losses and US $ 255 billion in lost labor income.
Ecosystems are also impacted by air pollution, particularly sulphur and
nitrogen emissions, and ground-level ozone as it affects their ability to
function and grow (Wang, 2021). Emissions of both sulphur dioxide and nitrogen
oxides deposit in water, on vegetation and on soils as ‘acid rain’, thereby
increasing their acidity with adverse effects on flora and fauna (De Vries et
al., 2014). Increased ground-level ozone also causes damage to cell membranes
on plants inhibiting key processes required for their growth and development
(Kim et al., 2020).
Ozone and PM also interact with radiation thus forcing climate change (Fiore
and Leibensperger, 2015; Denjean et al., 2016). PM warms by absorbing sunlight
(e.g., black carbon) or cools by scattering sunlight (e.g., sulfates) and
interacts with clouds; these radioactive and microphysical interactions can
induce changes in precipitation and regional circulation patterns.
All the above considerations clearly tell us that air quality is far from
being a ‘local’ problem: it affects human life in a myriad of different ways
and the same is for the whole Earth environment, a fact requiring effective,
global, pollution control/prevention measures. These latter are nowadays
possible thanks to the concomitant advent of supercomputers and of accurate
air pollution forecasting strategies. Observations, both from satellites and
from the ground, can provide a snapshot of the air quality, but have no real
predictive capability, at least on time horizons of the order of a few days.
Combining state-of-the-art numerical models of the atmosphere with satellite
and ground-station observations to provide daily forecasts of the composition
of the air worldwide, has become a key task to set up strategies of pollution
control and risk prevention. It is in this framework that the Air Quality
thematic area of CAMS, the Copernicus Atmosphere Monitoring Service, monitors
and forecasts European air quality long-range transport of pollutants (Marécal
et al., 2015). The combination of millions of daily observations and the
predictive power of Numerical Prediction Models (NPM) is a real strength of
CAMS. But this is not the only one. CAMS air quality forecasts indeed merge
the power of deterministic forecasts together with the power of statistical
forecasting. CAMS forecasts are indeed based on an ensemble of 9 state-of-the-
art numerical air quality models developed in Europe: CHIMERE from INERIS
(France), EMEP from MET Norway (Norway), EURAD-IM from Jülich IEK (Germany),
LOTOS-EUROS from KNMI and TNO (Netherlands), MATCH from SMHI (Sweden), MOCAGE
from Meteo-France (France), SILAM from FMI (Finland), DEHM from Aarhus
University (Denmark), and GEM-AQ from IEP-NRI (Poland). End-users have in this
way information about prediction uncertainty together with a single forecast
for different pollutants (corresponding, e.g., to the median of the empirical
distribution composed by the 9 members) for each forecast time horizon.
There is however a price to pay for obtaining worldwide air quality forecasts.
The necessarily coarse spatial resolution (of about 0.1o for CAMS) indeed
implies that small-scale features, including small-scale orography variations,
small-scale meteorological conditions and related dispersion regimes, local
pollutant emissions not captured by large-scale emission inventories, and many
others, are not explicitly accounted for in large-scale air-quality forecast
models. Using ‘raw’ CAMS forecasts to predict pollutant ground concentrations
at local scale may thus results in large forecast errors. The mean and
variance of the ensemble is expected to correlate with observations and the
actual model uncertainty, respectively, but a tendency to underestimate and to
be underdispersive with respect to the true observations and uncertainty of
the model is also expected. This is a reasonable expectation in analogy with
what happens for the 50-member-based meteorological Ensemble Prediction System
(EPS) of the European Centre for Medium-Range Weather Forecasts (ECMWF)
(Molteni et al., 1996; Montani et al., 2019). The lesson we learned from the
ECMWF EPS is that, used as is to predict, e.g., winds in regions of complex
orography, has a limited skill. On the contrary, accurate predictions can be
obtained via suitable statistical post processing methods of the raw
predictions which use past observed data to train the calibration algorithm
(Casciaro et al., 2022b).

Figure 1: The Liguria region (northwestern Italy) hosting the Regional
Monitoring Network used here to calibrate/validate the CAMS air quality
forecasts.

In this study, we explore the possibility of using CAMS air quality forecasts
to predict the 48-hour ahead evolution of PM10, PM2.5 (daily averages), O3,
NO2, and CO (hourly averages) in a network of stations located in Liguria
(northwestern Italy). The choice of the Ligurian region is motivated from the
fact that Liguria presents a large variety in its orographic shape (see Fig.
1), ranging from the alpine region (Alpi Marittime), with elevations up to
2000 meters, to flat areas as, e.g., in most zones close to the sea (Ligurian
Sea). Moreover, the mutual interaction between land and sea circulations makes
wind predictions, and predictions of weather variables in general, a very hard
problem. Because weather conditions heavily constrain air quality forecasts,
the same difficulties, probably even more severely, also affect pollutant
forecasts. As we will see, for the first 24-hour-ahead forecast horizon, the
added value of the raw CAMS forecasts with respect to persistence is very
limited (if there is one).
Here, we propose a statistical calibration of CAMS ensemble based on the use
of past pollutant measurements collected by the stations of the Ligurian
monitoring network (see Fig. 1) in the years 2019 and 2020. The calibration
strategy generalizes to air quality the one recently proposed by Casciaro et
al. (2022b) to calibrate the raw EPS forecasts for the wind speed and wind
power in complex orography regions for the green energy market.
As we will see, the calibration strategy remarkably increases the forecast
skill, thus making CAMS global forecasts appealing also for monitoring
activities at local scales.
Although our strategy has been applied to Liguria as an example, it can be
trivially extended to any other area in the world. For our calibration
strategy to be used operatively, the observed pollutant concentration on the
site of interest must be available in real/quasi-real time, other than as a
record of pairs of past observations/past forecasts.
The paper is organized as follow. In Sec. 2 a short introduction to the CAMS
multi-model ensemble system is presented; Sec. 3 presents the assessment of
raw forecasts of the CAMS multi-model system against observations collected by
the regional monitoring service of the Liguria region; Sec. 4 introduces our
calibration strategies, static and dynamic, while our results are shown in
Sec. 5. Concluding remarks and perspectives are drawn in the final section.

##  2 The CAMS air quality forecasts in a nutshell

The CAMS air quality models provide daily forecasts of the main atmospheric
pollutants concentrations on daily basis and with a look-ahead forecast
horizon of 4 days with a spatial resolution of 0.1 degrees and 7 vertical
levels from the Earth surface up to 5000 m. CAMS ensemble prediction system is
composed by 9 distinguishable air quality models, expression of state-of-the-
art knowledges in the field of multi-scale chemistry-transport models for
atmospheric composition. From the outputs of the 9 individual models, the
ensemble median can be calculated yielding forecast performances higher than
that of single models. The strength of CAMS ensemble lies in the fact that it
can provide an estimate of the forecast uncertainty via the spread between the
9 models.
The 9 models are distinguishable not only for the different details of the
physics/chemistry they describe but also for the fact that they rely on their
own data assimilation system. This latter consists of daily analyses of
pollutants near the surface by merging 1-day old observations with model
forecasts.
The meteorological fields needed to carry out the transport stage of
pollutants come from the 00:00 UTC operational ECMWF IFS forecast from the day
before, the boundary conditions for chemical species come from the CAMS IFS-
TM5 global production, while the emissions come from CAMS emission (for
anthropic emissions over Europe and for biomass burning).
CAMS forecasts are available for download from the CAMS Atmosphere Data Store
(ADS). The production time is 05:50 UTC (for 0-24 h forecasts) and 05:55 UTC
(for 24-48 h forecasts) with data availability guaranteed by 08:00 UTC, thus
being usable for operative needs.

Details on the multi-model forecast system are provided in Marécal et al.
(2015).

##  3 Assessing the raw CAMS air quality forecasts

Let us start our analysis by assessing the raw CAMS air quality forecasts
against measured data collected and validated by the monitoring service of the
Liguria region (see Fig. 1). The following pollutants have been considered in
the present study: PM10 (daily averages), PM2.5 (daily averages), O3 (hourly
averages), NO2 (hourly averages), and CO (hourly averages). Only 7 CAMS multi-
model ensemble members have been considered, GEM-AQ and DEHM becoming
available only from the last three months of 2019.

Table 1: Details on the observation stations grouped together by pollutant type. Pollutant |  | Number of
---
stations
| Data coverage
---
(%) 2019
| Data coverage
---
(%) 2020
CO | 13 | 56 | 73
NO2 | 34 | 74 | 63
O3 | 11 | 77 | 71
PM10 | 22 | 83 | 80
PM2.5 | 9 | 84 | 80

In Tab. 1 we report the number of ground stations for each of the measured
pollutants together with the total coverage for each of the two years
considered, 2019 and 2020 (i.e. the number of validated observations divided
by the maximum number of possible observations). The assessment of raw CAMS
forecasts is done in terms of standard error indices, both to assess the
accuracy of the mean (or median) of the ensemble and to quantify the accuracy
of the whole probability distribution of the multi-model ensemble forecasts.
The used indices are described in A.
Let us start by assessing the possible added value of CAMS forecasts against
the most economic and simple forecast one can do: the persistence based on
observations at a given time. When dealing with model evaluation this step is
strongly recommended. Here we selected the hour 9 UTC to build the
persistence-based forecasts. The results are reported in Fig. 2 where the
quality assessment of the CAMS ensemble mean against persistence is evaluated
in terms of the symmetrized skill score described in A. For the sake of
clarity, large positive values of the skill score mean that the CAMS forecasts
bring an added value with respect to persistence. The added value becomes more
and more important as the skill score approaches one. Negative values of the
skill score mean that, on average, persistence overcomes the raw CAMS
forecast. For CO, O3 and NO2 the look ahead time along the abscissa is
measured with respect to 9 UTC, the time used to define the persistence-based
forecasts. An abscissa of 1 thus means the look-ahead time corresponding to 10
UTC. For PM10 and PM2.5 the look ahead time is with respect to the analysis
time at 00 UTC.

Figure 2:  Mean over all stations of the skill score for both the NMAE index
and the correlation coefficient, ${\cal C}$, as a function of the look-ahead
time, for the mean of the raw CAMS ensemble forecasts The reference forecast
to calculate the skill score is based on persistence built from observations
at 9 UTC. Continuous lines: median of the skill scores computed by considering
all stations. Shaded areas represent the 50% confidence interval of the skill
score, accounting for the variability of the latter from station to station.

As one can clearly see from the figure, the raw CAMS multi-model ensemble
forecasts have a little added value with respect to persistence. The
conclusion holds more or less for all pollutants here considered, more
severely for CO and less dramatically for the PMs in the forecast interval
24-48 h. A complete panorama of the skill of the CAMS raw ensemble is
presented in B where all single raw members are assessed for the year 2020 by
averaging the single station indices, grouped together by pollutant type, over
all stations.

The main conclusion we can draw from this first step of our analysis is that
CAMS forecasts need to be calibrated using a training set of past observation-
forecast pairs. Before to discuss this important issue, let us try to
understand how CAMS multi-model ensemble performs with respect to the model
strategy, here referred to as S2020, used in Sartini et al. (2020) having
higher spatial resolution and accounting more detailed emission information
from the regional inventory emissions.

Table 2: A comparison between the hourly-averaged 24-h look ahead raw
forecasts from CAMS mean and from the high-resolution strategy reported in
Sartini et al. (2020). Error indices are calculated in both cases against one
year of observations of the Ligurian monitoring network.

Station name | NO2 | O3 | PM10
---|---|---|---
HH | NBI | HH | NBI | HH | NBI
Raw | S2020 | Raw | S2020 | Raw | S2020 | Raw | S2020 | Raw | S2020 | Raw | S2020
Mazzucca (IS) | 1.36 | 0.13 | -0.70 | 0.12 | - | - | - | - | 0.81 | - | -0.45 | -
Cengio (RS) | 0.59 | 0.51 | -0.31 | 0.26 | 0.45 | 0.41 | 0.48 | 0.15 | 0.48 | 0.56 | -0.09 | 0.35
Maggiolina (US) | 0.69 | 0.81 | -0.43 | 0.33 | 0.37 | 0.51 | 0.36 | 0.10 | 0.47 | 0.62 | -0.30 | 0.32
Quarto (US) | 0.53 | 0.83 | 0.20 | 0.43 | 0.28 | 0.34 | 0.16 | 0.07 | - | 0.47 | - | 0.31
Varaldo (US) | - | - | - | - | 0.15 | 0.47 | 0.05 | 0.08 | 0.61 | 0.75 | -0.22 | 0.34

To this aim, Tab. 2 reports the statistical error indices for a sub-set of
five control stations, classified as industrial (IS), rural (RS) and urban
(US) stations. We focused on such a sub-set so to have the possibility to
compare our results with those obtained by Sartini et al. (2020) using the
high-resolution numerical strategy, called LINEA, implemented and managed by
ARPA Liguria to forecast the concentration of photochemical pollutants. In the
table only daily means are considered as done in Sartini et al. (2020). The
added value brought by the concomitant action of high-resolution and more
detailed emission information is not clearly evident. Although a detailed
comparison among the two strategies is out of the scope of the present paper,
we can say that neither of the two strategies looks significantly better than
the other.
The question addressed here is the following: can a calibration strategy
properly set for ensemble prediction systems substantially increase the
forecast skill thereby overcoming persistence?
Answer this question is the main aim of the next sections of the present
paper. As a further aim we will also explore the possibility to obtain
calibrated forecasts significantly more accurate than those one can obtain
with costly high-resolution models which also account detailed emission
information made available at regional scale. Both answers have clear impacts
on the selection of optimal environmental monitoring strategies where the
optimality has also to account for the economic cost of the strategy.

##  4 A dynamic calibration strategy for multi-model ensemble

The calibration strategy used in the present activity has its roots in the
well-established Ensemble Model Output Statistics (EMOS) through which a raw
ensemble forecast can be accurately transformed into a predictive probability
density function, with a simultaneous correction of biases and dispersion
errors (Gneiting et al., 2005; Thorarinsdottir and Gneiting, 2010). If on the
one hand the technique is widely used to calibrate ensemble meteorological
models with remarkable improvements with respect to raw forecasts, the same
technique is still in its infancy in the realm of air quality modelling. The
reason is probably due to the fact that if on the side of meteorological
models several high-quality ensemble forecast systems are operative in many
leading weather centres, the same situation does not apply for operative
ensemble systems focused on air quality forecasts.
The technique we exploit here is largely inspired by a recent work of ours
where it was presented and applied to calibrate wind speed at the SYNOP
stations covering the Italian territory, with excellent results (Casciaro et
al., 2022a). The strategy is intimately dynamics in that it uses real
time/quasi real time observed data at station to reduce the real-time model
error and generalizes the static approach proposed in a previous paper by
Casciaro et al. (2022b) where observed data only enter in a training set of
past forecast-observation pair at stations.
To summarize the core of our strategy, called EMOS+4r, let us consider
$X_{i}^{j}$ the i-th ensemble member forecast of K members on the j-th model
grid-point ($j=1,\cdots,4$ spans over the nearest model grid-points to the
ground station) and $S^{2\>j}$ the members variance on the j-th model grid-
point. We also denote by $Z_{1},\cdots,Z_{q}$ the $q$ categorical variables
(here, the 2-m temperature, the persistence built from the last available
pollutant observation, the CAMS ensemble mean for the predictant pollutant,
the hour of the day, the surface wind gust, and the 10-m wind speed) expected
to be useful to disentangle the forecast error. Meteorological conditioning
variables are obtained here from the WRF numerical weather prediction model
(Skamarock et al., 2005) we have run for the years 2019 and 2020 in the model
configuration detailed in Cassola et al. (2015).
The EMOS+4r strategy defines a predictive probability distribution, a gamma
distribution in our strategy, with mean and variance given by

$\mu=a(Z_{1},\cdots,Z_{q})+\sum_{i=1,j=1}^{M,4}b_{ij}(Z_{1},\cdots,Z_{q})X_{i}^{j}$
(1)
$\sigma^{2}=c(Z_{1},\cdots,Z_{q})+\sum_{j=1}^{4}d_{j}(Z_{1},\cdots,Z_{q})S^{2\>j}$
(2)

where $j$ spans over the 4 model grid points around the station.
The free parameters $a$, $b$, $c$, and $d$, actually non-parametric functions
of the conditioning variables, must be best-fitted for each combination of
classes’ levels, via a training set (the year 2019 in the present study), by
minimizing the ‘distance’ between observed data and the associated probability
of occurrence. In plain terms, the quantity to be minimized is the so-called
Continuous Ranked Probability Score (CRPS) (Hersbach, 2000; Gneiting et al.,
2005).
For a gamma distribution, a closed form for the CRPS has been obtained by
Scheuerer and Möller (2015) making the minimization procedure easy and fast.
For an observation-forecast pair ($Y,\mathbf{X}$) it reads:

$crps=Y\left[2P\left(k,\frac{Y}{\theta}\right)-1\right]-k\theta\left[2P\left(k+1,\frac{Y}{\theta}\right)-1\right]-\frac{\theta}{\beta\left(\frac{1}{2},k\right)}$
(3)

with $Y$ being the observation, $P$ the incomplete gamma function (Abramowitz
and Stegun, 1948), and $\beta$ the beta function. The forecast vector
$\mathbf{X}=(X_{1},\cdots,X_{M})$ comes into the expression (3) via the
parameters $k$ and $\theta$. The quantity to be minimized in a training set
where both observations and forecasts are available is:

$CRPS=\frac{1}{N}\sum_{i=1}^{N}crps(\mathbf{X}_{i},Y_{i})$ (4)

with $i$ denoting the i-th observation-forecast pair and $N$ is the total
number of pairs in the training set.
The CRPS combines calibration and informativeness in one index, thus allowing
the evaluation of predictive performance that is based on the paradigm of
maximizing the sharpness of the predictive distributions subject to
calibration (Gneiting et al., 2007).
The EMOS+4r static calibration is completed by applying downstream to the
aforementioned steps a further EMOS with a 40-day rolling training set without
conditioning on $\mathbf{Z}$.
How to add a dynamic character to the present strategy is described in detail
by Casciaro et al. (2022a). Here we only remember that the basic idea of the
strategy is to insert two different persistences as new additional predictors
in (1) and (2). One persistence is built from station observation at hour 9
UTC (the choice for this specific hour can be modified as desired) while the
second type of persistence accounts for the possible occurrence of diurnal
cycles. For a leading time h, it is thus defined as the observed data at the
hour h-24. Using persistences as new predictors allows one to synchronize
forecasts to the nearest observed data with benefits visible in the first few
forecast hours as we will detail in the following.
We call our dynamic calibration D-EMOS+4r which exactly follows the steps
reported by Casciaro et al. (2022a).

##  5 Results

###  5.1 Choice of conditioning variables

The following set of conditioning variables has been considered in our study
for the predictands O3, NO2 and CO: the 2-m temperature (T2m), the persistence
using the predictand observations at 9 UTC (Pers), the mean of the CAMS
predictand forecasts (Ens mean), the hour of the day (h), the surface wind
gust (wg), and the 10-m wind speed (ws10).
For each station, the most effective three conditioning variables (i.e. M=3)
for which the minimum CRPS is obtained are determined via a cross validation.
This latter has been carried out on the training set corresponding to the year
2019, with the year 2020 serving as test set.

Figure 3: Shown is the frequency of occurrence of the conditioning variables
for CO, NO2 and O3. The conditioning variables we have considered are: the 2-m
temperature (T2m), the persistence using the predictand observation at 9 UTC
(Pers), the mean of the CAMS predictand forecasts (Ens mean), the hour of the
day (h), the surface wind gust (wg), and the 10-m wind speed (ws10).

Figure 3 shows the frequency with which each conditioning variable is
selected. A given percentage means the number of stations where a certain
conditioning variable has been selected divided by the number of stations
considered for the specific pollutant. From the figure it appears that the
most important conditioning variable is the 2-m temperature signaling a role
of the thermal structure of the atmospheric boundary layer in disentangling
the forecast error. The same seems to be confirmed by the important role of
the hour of the day as conditioning variable.
For the PMs, the importance of conditioning variables turned out to be
negligible. For this reason, the calibration strategy for the PMs will proceed
by assuming M=0 in the algorithm.

###  5.2 Assessing the added value of our calibration

Having selected the most effective three conditioning variables for O3, NO2
and CO, we are now ready to assess the skill of our calibrated forecasts using
the year 2020 as test set. Let us start by considering the simplest static
EMOS strategy where neither the conditioning variables nor the rolling
training set are considered, and the forecast is for the model grid point
closest to the station. Let us call this baseline calibration EMOS0. The skill
of the EMOS0-based calibration is shown in Fig. 4 using 9 UTC persistence as a
reference.

Figure 4: As in Fig. 2 but now the skill of the EMOS0-based calibration is
assessed against persistence built from the 9 UTC observation.

Although improvements with respect to the raw forecasts shown in Fig. 2 are
clearly detectable, persistence is still the best forecast for CO, for almost
all look-ahead times, and for the shortest look-ahead times for all other
pollutants.
These remarks suggest the use of more elaborated calibrations. Let us start
from our static EMOS+4r calibration the added value of which, against
persistence, is shown in Fig. 5.

Figure 5:  As in Fig. 2 but now the skill of our static calibration EMOS+4r is
assessed against persistence.

A further improvement is clearly observed downstream of the EMOS+4r
calibration even if problems still remain. In particular, for CO persistence
remains the best choice for all forecast horizons even if it is losing its
lead. In all cases persistence continues to overcome the calibrated forecasts
for the shortest forecast horizons.
This last remark strongly suggests the use of real-time data to correct the
model error in its infancy. This is done here in terms of our dynamic
calibration D-EMOS+4r the skill of which is shown in Fig. 6 using 9 UTC
persistence as a reference.

Figure 6:  As in Fig. 2 but now the skill of our dynamic calibration D-EMOS+4r
is assessed against 9 UTC persistence.

A net change of paradigm brought by our dynamic calibration D-EMOS+4r is
detectable. Persistence is no longer the best choice for CO and the calibrated
forecast increases its lead over persistence for all pollutants considered.
To better isolate the added value of our dynamic calibration D-EMOS+4r against
our static EMOS+4r calibration, in Fig. 7 we have reported the skill of the
D-EMOS+4r-based forecasts using the EMOS+4r-based forecasts as a reference.

Figure 7:  As in Fig. 2 but now the skill of our dynamic calibration D-EMOS+4r
is assessed against static calibration EMOS+4r.

The improvement of D-EMOS+4r relative to EMOS+4r is evident for all look-ahead
horizons, being more pronounced, as expected, for the shortest forecast
horizons.
Our best calibration D-EMOS+4r can be finally assessed with respect to the raw
forecasts in Fig. 8.

Figure 8:  As in Fig. 2 but now the skill of our dynamic calibration D-EMOS+4r
is assessed against raw forecast (mean of CAMS members).

The enhanced quality of the calibrated forecasts with respect to the raw
forecasts is noteworthy. A similar conclusion holds for the assessment of
calibration and informativeness encoded in the indices CRPS and $\Delta$ (see
A), as one can see from Fig. 9.

Figure 9:  As in Fig. 2 but now the skill of D-EMOS+4r is assessed in terms of
indices $\Delta$ and CRPS using the raw forecasts (CAMS members mean) as
reference.

From all figures shown in this section we can answer the question we addressed
in Sec. 3 on whether a calibration strategy properly set for ensemble
prediction systems can substantially increase the skill of air quality multi-
model forecasts. The answer is definitely positive and a key role played by
real-time observations clearly emerge.

We are now ready to answer the last question addressed in Sec. 3. It can be
rephrased here as follow. Is our best calibration competitive, or even better,
when compared against forecasts produced by raw high-resolution air quality
models which also account for detailed emission information made available at
regional scale?

Table 3: As in Tab. 2 but now the mean of the raw CAMS members is replaced by
our dynamic calibration D-EMOS+4r (D-EMOS in short).

Station name | NO2 | O3 | PM10
---|---|---|---
HH | NBI | HH | NBI | HH | NBI
D-EMOS | S2020 | D-EMOS | S2020 | D-EMOS | S2020 | D-EMOS | S2020 | D-EMOS | S2020 | D-EMOS | S2020
Mazzucca (IS) | 0.25 | 0.13 | 0.04 | 0.12 | - | - | - | - | 0.28 | - | 0.02 | -
Cengio (RS) | 0.32 | 0.51 | 0.00 | 0.26 | 0.18 | 0.41 | 0.02 | 0.15 | 0.39 | 0.56 | -0.01 | 0.35
Maggiolina (US) | 0.23 | 0.81 | 0.03 | 0.33 | 0.17 | 0.51 | 0.04 | 0.10 | 0.26 | 0.62 | 0.03 | 0.32
Quarto (US) | 0.39 | 0.83 | 0.09 | 0.43 | 0.17 | 0.34 | 0.03 | 0.07 | - | 0.47 | - | 0.31
Varaldo (US) | - | - | - | - | 0.13 | 0.47 | -0.01 | 0.08 | 0.35 | 0.75 | -0.03 | 0.34

Table 3 gives us a clear answer. Shown in the table are all error indices we
have computed for the five control stations already considered in Tab. 2. The
table also includes the error indices for the same stations associated to the
forecast model used by Sartini et al. (2020). The conclusion is that the
potential benefits brought by the combined use of a high-resolution model and
a more detailed inventory information can be significantly exceeded upon
implementing a proper dynamic calibration of the coarser model strategy.

We conclude this section by summarizing the skill of the calibrated CAMS
ensemble in terms of all considered indices. These latter have been calculated
as a mean of the indices over all stations. The results grouped together by
pollutant type are reported in Tab. 4.

Table 4: Summary of all relevant statistical error indices for the mean of the 9 CAMS members on the year 2020. Indices have been obtained by averaging all single index over all stations. Daily averages are considered for PM10 and PM2.5 while hourly averages for O3 and NO2 and CO. |  | Raw ens.
---
mean
Persistence |  | Static
---
calibration
| Dynamic
---
calibration
NMAE | 0.71 | 0.38 | 0.38 | 0.25 | CO
$\cal{C}$ | 0.34 | 0.63 | 0.41 | 0.75
HH | 1.57 | 0.48 | 0.51 | 0.34
NBI | -0.69 | 0.20 | -0.02 | -0.02
$\Delta$ | 1.61 | - | 0.32 | 0.19
CRPS | 0.39 | - | 0.15 | 0.10
NMAE | 0.27 | 0.28 | 0.21 | 0.18 | O3
$\cal{C}$ | 0.73 | 0.62 | 0.75 | 0.81
HH | 0.30 | 0.38 | 0.25 | 0.22
NBI | 0.15 | -0.13 | 0.01 | 0.01
$\Delta$ | 0.48 | - | 0.15 | 0.14
CRPS | 13.2 | - | 9.5 | 8.4
NMAE | 0.67 | 0.63 | 0.42 | 0.36 | NO2
$\cal{C}$ | 0.41 | 0.41 | 0.60 | 0.70
HH | 1.17 | 0.69 | 0.51 | 0.45
NBI | -0.39 | 0.29 | 0.04 | 0.03
$\Delta$ | 0.98 | - | 0.13 | 0.10
CRPS | 11.3 | - | 5.6 | 4.8
NMAE | 0.41 | 0.34 | 0.34 | 0.28 | PM2.5
$\cal{C}$ | 0.53 | 0.63 | 0.55 | 0.72
HH | 0.51 | 0.47 | 0.48 | 0.39
NBI | -0.14 | 0.01 | 0.00 | 0.01
$\Delta$ | 0.51 | - | 0.20 | 0.19
CRPS | 3.9 | - | 2.8 | 2.4
NMAE | 0.42 | 0.28 | 0.27 | 0.22 | PM10
$\cal{C}$ | 0.51 | 0.58 | 0.54 | 0.69
HH | 0.63 | 0.39 | 0.36 | 0.31
NBI | -0.33 | 0.00 | 0.02 | 0.02
$\Delta$ | 0.85 | - | 0.20 | 0.19
CRPS | 6.8 | - | 3.7 | 3.1

##  6 Conclusions and perspectives

The CAMS air quality model has been assessed for PM10, PM2.5, O3, NO2, and CO
against observations collected by the Regional Monitoring Network of the
Liguria region (northwestern Italy) in the years 2019 and 2020. Liguria has a
complex orography and is characterized by mutual interactions between land and
sea circulations making accurate weather predictions a difficult task. The
same holds true for air quality forecasts because of the key role of weather
in conditioning air quality. The selected observation network is thus a
particularly selective benchmark to assess air quality forecast strategies.
The main conclusion of this first step of our analysis is that in its raw form
the CAMS predictive strength is very limited, with the forecast built from
persistence having higher skills for all pollutants considered.
Motivated by these evidences we have asked the question on whether a suitable
calibration strategy acting as a post-processing of the raw CAMS forecasts may
bring out the intrinsic strength of CAMS multi-model ensemble strategy. To
answer this question we have applied a dynamic calibration strategy, named
D-EMOS+4r, from which a positive answer to our question definitely came. Our
results also highlight the key role of pollutant real-time observations to be
ingested in the calibration strategy. They bring a clear added value
especially in the shorter look-ahead forecast hours.
These results also point out the key value of observations, including their
management to be used in real time: they serve not only to have detailed
information on the current air quality state but also to predict with accuracy
the future scenario with positive outcomes for the reliability of pollution
early warning systems.

A second question was also answered in relation to the possible added value of
higher-resolution forecast strategies ingesting detailed information from the
inventory emission. When comparing the raw CAMS forecasts against the higher-
resolution strategy by Sartini et al. (2020) the aforementioned added value
hardly emerges. The situation of substantial equivalence is however totally
washed out when our dynamic calibration D-EMOS+4r is applied to the raw CAMS
forecasts, thus bringing out the intrinsic strength of CAMS multi-model
ensemble strategy. Left to future investigations, it will be interesting to
see whether our strategy applied to higher-resolution forecasts brings
additional advantages.
Our new dynamic calibration gave birth to the operative air-quality prediction
system for the Liguria region (DICCA, 2022).

##  7 Acknowledgments

G.C. has been funded by the Italian bank foundation “Fondazione Carige”. A.M.
acknowledges the funding from the Interreg Italia-Francia Marittimo SICOMAR+
project (grant number D36C17000120006) and from the Compagnia di San Paolo
(Project MINIERA No. I34I20000380007). Discussions with Francesco Ferrari and
Ludovica Sartini are warmly acknowledged. Observational data have been
provided by ARPAL the activity of which is acknowledged.

##  Appendix A Statistical indices

The following error indices are considered in the present study: the
normalized mean absolute error (NMAE), the normalized bias (NBI), the
symmetrically normalized root mean square error (HH), the correlation
coefficient (${\cal C}$), and the so-called reliability index ($\Delta$)
proposed by Delle Monache et al. (2006). The advantage of HH with respect to
the widely used non-symmetrically normalized root mean square error, RMSE, is
extensively discussed by Mentaschi et al. (2013).

By denoting with $Y_{n}$ the n-th observation, $X_{n}$ the corresponding n-th
forecast (here corresponding to the mean of the CAMS ensemble) and $N$ the
number of observation-forecast pairs in a given test set (here the year 2020),
the above mentioned statistical indices are defined as:

$NMAE=\frac{\sum_{n=1}^{N}\left|X_{n}-Y_{n}\right|}{\sum_{n=1}^{N}Y_{n}}$ (5)
$NBI=\frac{\sum_{n=1}^{N}\left(Y_{n}-X_{n}\right)}{\sum_{n=1}^{N}Y_{n}}$ (6)
$HH=\left(\frac{\sum_{n=1}^{N}\left(Y_{n}-X_{n}\right)^{2}}{\sum_{n=1}^{N}Y_{n}X_{n}}\right)^{1/2}$
(7)
${\cal{C}}=\frac{\sum_{n=1}^{N}(X_{n}-\overline{X})(Y_{n}-\overline{Y})}{N\sigma_{X}\sigma_{Y}}$
(8)

with:

$\sigma_{X}=\sqrt{\frac{\sum_{n=1}^{N}(X_{n}-\overline{X})^{2}}{N}}$ (9)
$\sigma_{Y}=\sqrt{\frac{\sum_{n=1}^{N}(Y_{n}-\overline{Y})^{2}}{N}}$ (10)

where $\overline{X}$ and $\overline{Y}$ are the mean values of $X$ and $Y$.
The goal of the probabilistic forecast, according to Gneiting et al. (2007),
is to maximize the sharpness of the predictive distribution subject to
calibration. Anderson (1996) and Hamill and Colucci (1997) proposed the use of
verification rank (VR) histograms to assess the calibration of ensemble
forecasts. VR histograms show the distribution of the ranks when the ranks of
the observations are pooled within the ordered ensemble forecasts. In a
calibrated ensemble, the observations and ensemble predictions should be
interchangeable, resulting in a uniform VR histogram. The continuous analogue
of the VR histogram is the probability integral transform (PIT) histogram
(Dawid, 1984; Diebold et al., 1997; Gneiting et al., 2007). The PIT value is
determined by the value of the predictive cumulative distribution function at
the verifying observation. For calibrated forecasts, the empirical cumulative
distribution function of PIT values should converge to the uniform
distribution.
The reliability index $\Delta$ was proposed by Delle Monache et al. (2006) to
quantify the deviation of VR histograms from uniformity. To quantify the
deviation from uniformity in the PIT histograms, we use here the following
definition of $\Delta$:

$\Delta=\sum_{i=1}^{m}\left|f_{i}-\frac{1}{m}\right|$ (11)

where $m$ is the histogram number of classes, each with a relative frequency
of $1/m$, and $f_{i}$ is the observed relative frequency in class $i$.
This index ranges from 0 to $+\infty$, with the bound zero corresponding to
optimal forecast.

Having defined the statistical error indices, a symmetric variant of the
standard Skill Score index (Wilks, 2011) is used here to make the comparison
between different forecasts (raw or calibrated) as quantitative as possible.
Our new definition is used to assess the performance of a given forecast by
comparing its associated statistical error index against the one corresponding
to a reference forecast. More quantitatively, the symmetrized skill score, SS,
is defined here as

$SS=\frac{A-A_{ref}}{A_{opt}-A_{ref}}\qquad\mbox{when}\qquad\frac{A-A_{ref}}{A_{opt}-A_{ref}}\geq
0$ $SS=-\frac{A_{ref}-A}{A_{opt}-A}\qquad\mbox{elsewhere}$

where $A$ is the value of a suitable error index associated to forecast,
$A_{ref}$ is the same as $A$ but relative to a reference forecast. Finally,
$A_{opt}$ refers to the optimal index value here supposed to be 0 (e.g. as for
NMAE) or 1 (e.g. as for the correlation coefficient). All these error indices
are supposed here to be positive (e.g. as it is for the NMAE and for the
correlation coefficient apart situations of anticorrelated regimes however not
of interest here).

A perfect calibration yields SS=1, corresponding to the upper bound of SS.
Values of SS smaller than one (including negative values) indicate that the
forecast is less accurate than the reference one.

##  Appendix B Assessment of the skill of the single CAMS members for the year
2020

For the sake of comparison, we assess here the skill of the 9 CAMS ensemble
members. Error indices have been calculated for the year 2020 for PM10, PM2.5
(daily averages), O3, NO2, and CO (hourly averages) against observations from
the Ligurian monitoring network. CAMS median and mean are also assessed
together with the 9 members. While for the majority of indices it is hard to
identify a model clearly overcoming the others, for the correlation
coefficient CAMS mean and median show the highest skill, thus confirming the
power of the ensemble approach.

Table 5: Statistical error indices for the 9 raw CAMS members: CHIMERE (CH)
from INERIS (France), EMEP (EM) from MET Norway (Norway), EURAD-IM (EU) from
Jülich IEK (Germany), LOTOS-EUROS (LO) from KNMI and TNO (Netherlands), MATCH
(MA) from SMHI (Sweden), MOCAGE (MO) from Meteo-France (France), SILAM (S)
from FMI (Finland), DEHM (DE) from Aarhus University (Denmark), and GEM-AQ
(GE) from IEP-NRI (Poland). The ensemble mean and median are, Mean and Median,
respectively.

| EM | SI | MO | EU | LO | MA | CH | GE | DE | Median | Mean |
---|---|---|---|---|---|---|---|---|---|---|---|---
NMAE | 0.69 | 0.68 | 0.70 | 0.73 | 0.71 | 0.73 | 0.71 | 0.69 | 0.73 | 0.71 | 0.71 | CO
HH | 1.5 | 1.4 | 1.6 | 1.7 | 1.6 | 1.7 | 1.6 | 1.5 | 1.7 | 1.6 | 1.6
NBI | -0.66 | -0.66 | -0.68 | -0.72 | -0.69 | -0.72 | -0.69 | -0.66 | -0.71 | -0.70 | -0.69
$\cal{C}$ | 0.17 | 0.31 | 0.29 | 0.20 | 0.24 | 0.19 | 0.22 | 0.21 | 0.23 | 0.26 | 0.27
NMAE | 0.29 | 0.31 | 0.41 | 0.36 | 0.35 | 0.36 | 0.38 | 0.37 | 0.36 | 0.30 | 0.29 | O3
HH | 0.33 | 0.39 | 0.42 | 0.42 | 0.37 | 0.37 | 0.39 | 0.41 | 0.38 | 0.32 | 0.31
NBI | 0.12 | -0.09 | 0.31 | 0.11 | 0.22 | 0.28 | 0.28 | 0.21 | 0.25 | 0.20 | 0.19
$\cal{C}$ | 0.66 | 0.63 | 0.58 | 0.41 | 0.62 | 0.65 | 0.60 | 0.58 | 0.57 | 0.71 | 0.73
NMAE | 0.68 | 0.69 | 0.78 | 0.71 | 0.73 | 0.74 | 0.76 | 0.75 | - | 0.68 | 0.67 | NO2
HH | 1.3 | 0.99 | 1.1 | 1.4 | 1.3 | 1.7 | 1.3 | 1.3 | - | 1.3 | 1.2
NBI | -0.43 | -0.17 | -0.13 | -0.47 | -0.43 | -0.63 | -0.34 | -0.34 | - | -0.43 | -0.37
$\cal{C}$ | 0.33 | 0.38 | 0.38 | 0.31 | 0.28 | 0.30 | 0.16 | 0.23 | - | 0.37 | 0.38
NMAE | 0.58 | 0.52 | 0.47 | 0.48 | 0.55 | 0.44 | 0.46 | 0.57 | 0.46 | 0.41 | 0.40 | PM2.5
HH | 0.72 | 0.64 | 0.84 | 0.67 | 0.93 | 0.69 | 0.58 | 0.72 | 0.62 | 0.61 | 0.57
NBI | 0.11 | 0.16 | -0.39 | -0.07 | -0.41 | -0.24 | 0.05 | 0.21 | -0.11 | -0.13 | -0.08
$\cal{C}$ | 0.43 | 0.50 | 0.51 | 0.45 | 0.40 | 0.50 | 0.49 | 0.38 | 0.55 | 0.53 | 0.55
NMAE | 0.46 | 0.39 | 0.61 | 0.48 | 0.54 | 0.48 | 0.44 | 0.48 | 0.42 | 0.44 | 0.41 | PM10
HH | 0.65 | 0.55 | 1.14 | 0.77 | 0.83 | 0.75 | 0.65 | 0.74 | 0.61 | 0.59 | 0.63
NBI | -0.22 | -0.17 | -0.60 | -0.41 | -0.32 | -0.41 | -0.29 | -0.29 | -0.31 | -0.37 | -0.34
$\cal{C}$ | 0.44 | 0.52 | 0.49 | 0.47 | 0.41 | 0.50 | 0.43 | 0.32 | 0.56 | 0.54 | 0.55

## References

  * Abramowitz and Stegun (1948) Abramowitz, M., Stegun, I. A., 1948. Handbook of mathematical functions with formulas, graphs, and mathematical tables. Vol. 55. US Government printing office.
  * Anderson (1996) Anderson, J. L., 1996. A method for producing and evaluating probabilistic forecasts from ensemble model integrations. Journal of Climate 9 (7), 1518–1530.
  * Casciaro et al. (2022a) Casciaro, G., Ferrari, F., Lagomarsino Oneto, D., Lira Loarca, A., Mazzino, A., 2022a. Increasing the skill of short-term wind speed ensemble forecasts combining forecasts and observations via a new dynamic calibration.
  * Casciaro et al. (2022b) Casciaro, G., Ferrari, F., Mazzino, A., 2022b. Novel strategies of Ensemble Model Output Statistics (EMOS) for calibrating wind speed/power forecasts.
  * Cassola et al. (2015) Cassola, F., Ferrari, F., Mazzino, A., 2015. Numerical simulations of Mediterranean heavy precipitation events with the WRF model: A verification exercise using different approaches. Atmospheric Research 164-165, 210–225.
URL https://www.sciencedirect.com/science/article/pii/S0169809515001696

  * Dawid (1984) Dawid, A. P., 1984. Present position and potential developments: Some personal views statistical theory the prequential approach. Journal of the Royal Statistical Society: Series A (General) 147 (2), 278–290.
  * De Vries et al. (2014) De Vries, W., Dobbertin, M. H., Solberg, S., Van Dobben, H. F., Schaub, M., 2014\. Impacts of acid deposition, ozone exposure and weather conditions on forest ecosystems in Europe: an overview. Plant and Soil 380 (1), 1–45.
URL https://doi.org/10.1007/s11104-014-2056-2

  * Delle Monache et al. (2006) Delle Monache, L., Hacker, J. P., Zhou, Y., Deng, X., Stull, R. B., 2006\. Probabilistic aspects of meteorological and ozone regional ensemble forecasts. Journal of Geophysical Research: Atmospheres 111 (D24).
  * Denjean et al. (2016) Denjean, C., Cassola, F., Mazzino, A., Triquet, S., Chevaillier, S., Grand, N., Bourrianne, T., Momboisse, G., Sellegri, K., Schwarzenbock, A., Freney, E., Mallet, M., Formenti, P., 2016. Size distribution and optical properties of mineral dust aerosols transported in the western Mediterranean. Atmospheric Chemistry and Physics 16 (2), 1081–1104.
URL https://acp.copernicus.org/articles/16/1081/2016/

  * DICCA (2022) DICCA, 2022. The operative air-quality prediction system for the Liguria region. http://www3.dicca.unige.it/mazzino/Air-Quality-Forecasts/forecast.html.
  * Diebold et al. (1997) Diebold, F. X., Gunther, T. A., Tay, A., 1997. Evaluating density forecasts.
  * Fiore and Leibensperger (2015) Fiore, A., Leibensperger, E., 2015. Air quality and climate connections. Journal Air Waste Management Association 65 (6), 645–685.
  * Gneiting et al. (2007) Gneiting, T., Balabdaoui, F., Raftery, A. E., 2007. Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69 (2), 243–268.
  * Gneiting et al. (2005) Gneiting, T., Raftery, A. E., Westveld III, A. H., Goldman, T., 2005\. Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation. Monthly Weather Review 133 (5), 1098–1118.
  * Hamill and Colucci (1997) Hamill, T. M., Colucci, S. J., 1997. Verification of Eta–RSM short-range ensemble forecasts. Monthly Weather Review 125 (6), 1312–1327.
  * Hersbach (2000) Hersbach, H., 2000. Decomposition of the continuous ranked probability score for ensemble prediction systems. Weather and Forecasting 15 (5), 559–570.
  * Kim et al. (2020) Kim, J. J., Fan, R., Allison, L. K., Andrew, T. L., 2020. On-site identification of ozone damage in fruiting plants using vapor-deposited conducting polymer tattoos. Science Advances 6 (36), eabc3296.
URL https://www.science.org/doi/abs/10.1126/sciadv.abc3296

  * Marécal et al. (2015) Marécal, V., Peuch, V.-H., Andersson, C., Andersson, S., Arteta, J., Beekmann, M., Benedictow, A., Bergström, R., Bessagnet, B., Cansado, A., Chéroux, F., Colette, A., Coman, A., Curier, R. L., Denier van der Gon, H. A. C., Drouin, A., Elbern, H., Emili, E., Engelen, R. J., Eskes, H. J., Foret, G., Friese, E., Gauss, M., Giannaros, C., Guth, J., Joly, M., Jaumouillé, E., Josse, B., Kadygrov, N., Kaiser, J. W., Krajsek, K., Kuenen, J., Kumar, U., Liora, N., Lopez, E., Malherbe, L., Martinez, I., Melas, D., Meleux, F., Menut, L., Moinat, P., Morales, T., Parmentier, J., Piacentini, A., Plu, M., Poupkou, A., Queguiner, S., Robertson, L., Rouïl, L., Schaap, M., Segers, A., Sofiev, M., Tarasson, L., Thomas, M., Timmermans, R., Valdebenito, A., van Velthoven, P., van Versendaal, R., Vira, J., Ung, A., 2015. A regional air quality forecasting system over Europe: the MACC-II daily ensemble production. Geoscientific Model Development 8 (9), 2777–2813.
URL https://gmd.copernicus.org/articles/8/2777/2015/

  * Mentaschi et al. (2013) Mentaschi, L., Besio, G., Cassola, F., Mazzino, A., 2013. Problems in RMSE-based wave model validations. Ocean Modelling 72, 53–58.
  * Molteni et al. (1996) Molteni, F., Buizza, R., Palmer, T. N., Petroliagis, T., 1996. The ECMWF ensemble prediction system: Methodology and validation. Quarterly Journal of the Royal Meteorological Society 122 (529), 73–119.
  * Montani et al. (2019) Montani, A., Gascon, E., Hewson, T., Paccagnella, T., Pillosu, F., Pincini, G., 2019. Performance of ECMWF- and COSMO-based ensemble forecast systems for precipitation events over Italy. EMS Annual Meeting Abstracts Vol. 16, EMS2019-146.
  * Sartini et al. (2020) Sartini, L., Antonelli, M., Pisoni, E., Thunis, P., 2020. From emissions to source allocation: Synergies and trade-offs between top-down and bottom-up information. Atmospheric Environment: X 7, 100088.
URL https://www.sciencedirect.com/science/article/pii/S2590162120300289

  * Scheuerer and Möller (2015) Scheuerer, M., Möller, D., 2015. Probabilistic wind speed forecasting on a grid based on ensemble model output statistics. The Annals of Applied Statistics 9 (3), 1328–1349.
  * Shore (2019) Shore, S. A., 2019. The metabolic response to ozone. Frontiers in Immunology 10\.
URL https://www.frontiersin.org/article/10.3389/fimmu.2019.02890

  * Skamarock et al. (2005) Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Wang, W., Powers, J. G., 2005. A description of the advanced research WRF version 2\. Tech. rep., National Center For Atmospheric Research Boulder Co Mesoscale and Microscale.
  * Thorarinsdottir and Gneiting (2010) Thorarinsdottir, T. L., Gneiting, T., 2010. Probabilistic forecasts of wind speed: Ensemble model output statistics by using heteroscedastic censored regression. Journal of the Royal Statistical Society: Series A (Statistics in Society) 173 (2), 371–388.
  * Wang (2021) Wang, Z., 2021. Satellite-Observed Effects from Ozone Pollution and Climate Change on Growing-Season Vegetation Activity over China during 1982-2020. Atmosphere 12 (11).
URL https://www.mdpi.com/2073-4433/12/11/1390

  * Wilks (2011) Wilks, D. S., 2011. Statistical methods in the atmospheric sciences. Vol. 100\. Academic Press.
  * World Bank (2016) World Bank, 2016. The Cost of Air Pollution: Strengthening the Economic Case for Action.
URL https://openknowledge.worldbank.org/handle/10986/25013

  * Yin et al. (2017) Yin, P., He, G., Fan, M., Chiu, K. Y., Fan, M., Liu, C., Xue, A., Liu, T., Pan, Y., Mu, Q., Zhou, M., 2017. Particulate air pollution and mortality in 38 of China’s largest cities: time series analysis. BMJ 356.
URL https://www.bmj.com/content/356/bmj.j667


# Let’s reinvent subtraction

Peter G. Doyle

(Version 240 dated 25 March 2022
No Copyright††thanks:  The author(s) hereby waive all copyright and related or
neighboring rights to this work, and dedicate it to the public domain. This
applies worldwide. )

###### Abstract

Subtraction is a powerful technique for creating new bijections from old.
Let’s reinvent it! While we’re at it, let’s reinvent division as well.

##  1 Matchings

Write

$f:A\equiv B$

and say ‘$f$ matches $A$ with $B$’ to mean that we know a suitable bijection
$f$ from $A$ to $B$, together with its inverse

$f^{-1}:A\equiv B.$

Write

$A\equiv B$

and say ‘$A$ matches $B$’ to mean that we know (or know we _could_ know) some
$f:A\equiv B$. We have

$A\equiv A;\;\;\;A\equiv B\implies B\equiv A;\;\;\;A\equiv B\land B\equiv
C\implies A\equiv C.$

(We refrain from saying that $\equiv$ is an equivalence relation since it is
inherently time-dependent.)

We can add and multiply matchings:

$A\equiv B\land C\equiv D\implies A+C\equiv B+D,$

where $+$ denotes disjoint union, and

$A\equiv B\land C\equiv D\implies A\times C\equiv B\times D.$

##  2 Respectful subtraction

Before addressing subtraction in general, let’s begin with _respectful
subtraction_ , an important special case. It’s so simple that it hardly
deserves to be called subtraction.

######  Definition.

For $f:A+C\equiv B+D$, $g:C\equiv D$, write $g\ll f$ and say ‘$g$ respects
$f$’ if

$\forall x\in C\;(f(x)\in D\implies g(x)=f(x)).$

######  Proposition 1 (Respectful subtraction).

If

$f:A+C\equiv B+D;\;\;\;g:C\equiv D;\;\;\;g\ll f$

then

$f\setminus g:A\equiv B,$

where

$f\setminus g(x)=f(x){\ \tt{if}\ }f(x)\in B{\ \tt{else}\ }f(g^{-1}(f(x))).$

Moreover,

$f\setminus g\ll f$

and

$f\setminus({f\setminus g})=g.$

Proof: Without loss of generality, assume $f$ is the identity on $A+B=C+D$.
$g$ fixes $C\cap D$ and matches $C\setminus D$ to $D\setminus C$. $f\setminus
g$ fixes $A\cap B$ and, taking its cue from $g^{-1}$, matches $A\setminus
B=D\setminus C$ to $B\setminus A=C\setminus D$. $\quad\spadesuit$

##  3 Subtraction

######  Proposition 2 (Subtraction).

If

$f:A+C\equiv B+D;\;\;\;g:C\equiv D$

with $D$ finite, then

$f\setminus g:A\equiv B,$

where

$f\setminus g(x)=(y:=f(x);\;{\tt{while}\ }y\in D{\ \tt{do}\
}y:=f(g^{-1}(y));\;{\tt{return}\ }y).$

Proof: This goes way back—see [1]. $\quad\spadesuit$

When $g\ll f$ we’re back to respectful subtraction:

######  Proposition 3.

If $g\ll f$ then

$f\setminus g(x)=f(x){\ \tt{if}\ }f(x)\in B{\ \tt{else}\
}f(g^{-1}(f(x))).\quad\spadesuit$

The fact that $f\setminus g\ll f$ is general:

######  Proposition 4.

$f\setminus g\ll f$. $\quad\spadesuit$

Idempotence of subtraction characterizes respectfulness:

######  Proposition 5.

$f\setminus(f\setminus g)=g\;\iff\;g\ll f.\quad\spadesuit$

This gives us a closure operation:

######  Proposition 6.

$f\setminus(f\setminus(f\setminus g))=f\setminus g,$

so

$f\setminus(f\setminus(f\setminus(f\setminus g)))=f\setminus(f\setminus
g).\quad\spadesuit$

##  4 Inclusion-exclusion

Subtraction generalizes to inclusion-exclusion. Let $P$ be a poset with
$\\{q:q\leq p\\}$ finite for all $p$. Given a family of finite sets
$A_{p},p\in P$, let

$A_{\leq p}=\sum_{q\leq p}A_{q},$

etc.

######  Proposition 7 (Inclusion-exclusion).

$\forall p\;A_{\leq p}\equiv B_{\leq p}\;\;\implies\;\;\forall p\;A_{p}\equiv
B_{p}.$

Proof: By induction: Assuming

$\forall q<p\;A_{q}\equiv B_{q}$

(true if $p$ is minimal) we have

$A_{<p}\equiv B_{<p}.$

Subtract from

$A_{\leq p}\equiv B_{\leq p}$

to get

$A_{p}\equiv B_{p}.\quad\spadesuit$

To be more explicit, define the projection map

$\pi_{A}:\sum_{p}A_{p}\to P,\;\pi(x)=p\iff x\in A_{p}.$

######  Proposition 8.

If

$g_{p}:A_{\leq p}\equiv B_{\leq p}$

then

$f_{p}:A_{p}\equiv B_{p}$

where

$f_{p}(x)=F(p,x),$ $F(p,x)=(y:=g_{p}(x);\;q:=\pi_{B}(y);\;{\tt{return}\ }y{\
\tt{if}\ }q=p{\ \tt{else}\ }F(p,\bar{F}(q,y)));\;$
$\bar{F}(p,x)=(y:=g_{p}^{-1}(x);\;q:=\pi_{A}(y);\;{\tt{return}\ }y{\ \tt{if}\
}q=p{\ \tt{else}\ }\bar{F}(p,F(q,y))).$

Proof: This is what you get if you trace it through. $\quad\spadesuit$

##  5 Extreme division

Subtraction will get you a long way toward automating the process of
generating bijections. But sometimes you will want to divide, and that’s when
things can get scary. (Like the old Marchant mechanical calculators, which
would make a satisfying ‘chunk’ when you hit $+$ or $-$, but would make a
terrifying racket, with the carriage scurrying to and fro, when you hit the
${\overset{\mbox{\tiny Auto}}{\div}}$ key.)

The way to keep things under control is to make sure you are multiplying
polynomials (in one or many variables), with the polynomial you are dividing
by having a unique extreme monomial $\omega$ for some linear function on the
space of degrees (in other words a singleton monomial on the boundary of its
Newton polytope). In this case you can use ‘extreme division’, whereby you
recursively subtract the mapping based on multiplication by $\omega$.

Suppose

$F:A\times C\equiv B\times C$

and

$G:B\times C\equiv A\times C\;\;\;.$

(We may choose $G=F^{-1}$, but we don’t require this.) For $\omega\in C$,
define

$\operatorname{xdiv}((F,G),\omega)=(f,g)$

where $f,g$ are the partial functions $f$ on $A$ and $g$ on $B$ defined via
the mutual recursion equations

$f(x)=((y,z):=F((x,\omega));\,{\tt{while}\ }z\neq\omega{\ \tt{do}\
}((y,z):=F((g(y),z));\,{\tt{return}\ }y)$
$g(x)=((y,z):=G((x,\omega));\,{\tt{while}\ }z\neq\omega{\ \tt{do}\
}((y,z):=G((f(y),z));\,{\tt{return}\ }y)$

If both $f$ and $g$ are total, we say that the pair $(F,G)$ is _X-divisible
for $\omega$_, and that $\omega$ is an _extreme point_ for the pair $(F,G)$.
This terminology springs from the following proposition.

######  Proposition 9 (Extreme division).

If $A,B,C$ are multinomials, and $F,G$ match terms of $A\cdot C$ and $B\cdot
C$, then $(F,G)$ is X-divisible for any extreme monomial $\omega$ of $C$.

Proof: By induction. $\quad\spadesuit$

Extreme division is what Conway and Doyle [2] had in mind when they wrote,
‘There is more to division than repeated subtraction.’ This must be a great
truth, because its negation would appear to be at least as true.

##  6 Mode d’emploi

Contrary to what we may appear to be claiming in Proposition 9, the bijection
yielded by the $\mathrm{xdiv}$ algorithm may fail to be a ‘matching’, because
it may take too long to compute, or otherwise fail to qualify as ‘suitable’,
the admittedly slippery condition that we slipped in as part of the definition
of a matching. The same goes for bijections obtained by inclusion-exclusion.
This is why we have taken care to announce 7 and 9 as ‘Propositions’, yielding
bijections proposed for consideration as ‘matchings’. This is contrary to
mathematical custom, and wrong-headed, but useful nevertheless.

Now it will often happen that a slow quotient bijection can be speeded up
immensely by ‘memoizing’, meaning that values of $f$ and $g$ are automatically
saved so that they don’t get computed over and over. This in itself may make
the bijection ‘suitable’.

Better still is to be able to see ‘what the bijection is doing’, so that it
can be defined, and proven to be a suitable bijection, without reference to
its origin as a quotient.

Here’s a case in point. Everyone knows that

$\binom{n}{k}=\binom{n}{n-k},$

and if you ask why, they will either tell you to match a $k$-subset to the
complementary $(n-k)$ \- subset, or compute

$\binom{n}{k}=\frac{n!}{k!(n-k)!}=\frac{n!}{(n-k)!(n-(n-k))!}=\binom{n}{n-k}.$

Taking our motto to be ‘follow the algebra’, we recast the computation as

$\binom{n}{k}k!(n-k)!=n!=\binom{n}{n-k}(n-k)!k!=\binom{n}{n-k}k!(n-k)!,$

where every step is backed by a matching. Now divide. The xdiv algorithm
yields a very inefficient computation of a very simple matching (see the code
in the appendix):

    ([[0, 1], [2, 3, 4]], [[0, 1, 2], [3, 4]])
    ([[0, 2], [1, 3, 4]], [[0, 1, 3], [2, 4]])
    ([[0, 3], [1, 2, 4]], [[0, 2, 3], [1, 4]])
    ([[0, 4], [1, 2, 3]], [[1, 2, 3], [0, 4]])
    ([[1, 2], [0, 3, 4]], [[0, 1, 4], [2, 3]])
    ([[1, 3], [0, 2, 4]], [[0, 2, 4], [1, 3]])
    ([[1, 4], [0, 2, 3]], [[1, 2, 4], [0, 3]])
    ([[2, 3], [0, 1, 4]], [[0, 3, 4], [1, 2]])
    ([[2, 4], [0, 1, 3]], [[1, 3, 4], [0, 2]])
    ([[3, 4], [0, 1, 2]], [[2, 3, 4], [0, 1]])

This follow-the-algebra matching differs from only slightly from taking the
complementary set. It’s arguably better. Do you agree?

## Appendix

    """
    gauss 000 - use extreme division to match (n choose k) to (n choose n-k)
    """

    def xdiv(FG,omega):
      F,G=FG
      def f(x):
        y,z=F((x,omega))
        while z!=omega:
          y,z=F((g(y),z))
        return y
      def g(x):
        y,z=G((x,omega))
        while z!=omega:
          y,z=G((f(y),z))
        return y
      return [f,g]

    import itertools

    def makeintolist(a): return [x for x in a]
    def sublist(a,c): return [a[x] for x in c]
    def num(n): return [x for x in range(n)]

    def combinations(a,k):
      return [makeintolist(a) for a in itertools.combinations(a,k)]

    def complement(s,a): return [x for x in s if x not in a]
    def separate(s,a): return [a,complement(s,a)]

    """
    choose(n,k) lists the ways for separating n into pieces of size k and n-k
    """

    def choose(n,k):
      s=range(n)
      return [separate(s,a) for a in combinations(s,k)]

    """
    compress(l) replaces elements of l by their relative ranks
    """

    def compress(l):
      m=sorted(l)
      return [m.index(x) for x in l]

    """
    binom(sigma,k) maps the permutation sigma to choose(n,k) * k! * (n-k)!
    """

    def binom(sigma,k):
      a=sigma[:k]
      b=sigma[k:]
      return [[sorted(a),sorted(b)],[compress(a),compress(b)]]

    """
    monib is the inverse of binom
    """

    def monib(abcd):
      ((a,b),(c,d))=abcd
      return sublist(a,c)+sublist(b,d)

    """
    flipkl maps  choose(n,k) * k! * (n-k)! to  choose(n,k) * (n-k)! * k!

    We need this because we decided to cast division in terms of matchings
    between A*C and B*C rather than A*C and B*D
    """

    def flipkl(abcd):
      ((a,b),(c,d))=abcd
      return [[a,b],[d,c]]

    """
    Our F and G would be the same, except for the flipping.
    """

    def F(abcd):
      ((a,b),(c,d))=abcd
      k=len(a)
      l=len(b)
      return flipkl(binom(monib(abcd),l))

    def G(abcd):
      ((a,b),(c,d))=abcd
      k=len(a)
      l=len(b)
      return binom(monib(flipkl(abcd)),l)

    def omega(n,k): return [num(k),num(n-k)]

    def match(n,k): return xdiv([F,G],omega(n,k))

    """
    TESTING
    """

    def column(l):
      print(*l,sep=’\n’)

    n=5
    k=2

    A=choose(n,k)
    (f,g)=match(n,k)
    B=[f(ab) for ab in A]
    column(zip(A,B))
    C=[g(ab) for ab in B]
    print(A==C)

## References

  * [1] Peter G. Doyle.  A category for bijective combinatorics, 2019, arXiv:1907.09015 [math.CO].  https://arxiv.org/abs/1907.09015.
  * [2] Peter G. Doyle and John Horton Conway.  Division by three, 1994, arXiv:math/0605779 [math.LO].  http://arxiv.org/abs/math/0605779.