# Long-term trends of light pollution assessed from SQM measurements and an
empirical atmospheric model††thanks: Generated using Copernicus Climate Change
Service (C3S) Information and Atmosphere Monitoring Service (CAMS)
Information, 2019.

Johannes Puschnig,1 Stefan Wallner,2,3 Axel Schwope,4 Magnus Näslund5
1Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71,
D-53121 Bonn, Germany
2ICA, Slovak Academy of Sciences, Dubravska cesta 9, 84503 Bratislava, Slovak
Republic
3Universität Wien, Institut für Astrophysik, Türkenschanzstraße 17, A-1180
Wien, Austria
4Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482
Potsdam, Germany
5Department of Astronomy, Stockholm University, AlbaNova University Centre,
SE-10691 Stockholm, Sweden
E-mail<EMAIL_ADDRESS>
(Accepted 2022 October 15. Received 2022 October 11; in original form 2022
August 31)

###### Abstract

We present long-term (4-10 years) trends of light pollution observed at 26
locations, covering rural, intermediate and urban sites, including the three
major European metropolitan areas of Stockholm, Berlin and Vienna. Our
analysis is based on i) night sky brightness (NSB) measurements obtained with
Sky Quality Meters (SQMs) and ii) a rich set of atmospheric data products. We
describe the SQM data reduction routine in which we filter for moon- and
clear-sky data and correct for the SQM “aging” effect using an updated version
of the twilight method of Puschnig et al. (2021). Our clear-sky, aging-
corrected data reveals short- and long-term (seasonal) variations due to
atmospheric changes. To assess long-term anthropogenic NSB trends, we
establish an empirical atmospheric model via multi-variate penalized linear
regression. Our modeling approach allows to quantitatively investigate the
importance of different atmospheric parameters, revealing that surface albedo
and vegetation have by far the largest impact on zenithal NSB. Additionally,
the NSB is sensitive to black carbon and organic matter aerosols at urban and
rural sites respectively. Snow depth was found to be important for some sites,
while the total column of ozone leaves impact on some rural places. The
average increase in light pollution at our 11 rural sites is 1.7% per year. At
our nine urban sites we measure an increase of 1.8% per year and for the
remaining six intermediate sites we find an average increase of 3.7% per year.
These numbers correspond to doubling times of 41, 39 and 19 years. We estimate
that our method is capable of detecting trend slopes shallower/steeper than
$\pm$1.5% per year.

###### keywords:

light pollution – atmospheric effects – techniques: photometric

††pubyear: 2022††pagerange: Long-term trends of light pollution assessed from
SQM measurements and an empirical atmospheric model††thanks: Generated using
Copernicus Climate Change Service (C3S) Information and Atmosphere Monitoring
Service (CAMS) Information, 2019.–18

##  1 Introduction

During the last decade, an ever increasing number of studies found evidence
that artificial light at night (ALAN) leads to negative consequences – not
only for astronomy – but also for ecosystems (e.g. Longcore & Rich, 2004;
Perkin et al., 2011), biodiversity (e.g. Hölker et al., 2010), animals (e.g.
Eisenbeis, 2006; Perkin et al., 2014; Mathews et al., 2015; Owens et al.,
2020; Parkinson et al., 2020) and human beings (e.g. Chepesiuk, 2009; Haim &
Portnov, 2013; Cho et al., 2015; Garcia-Saenz et al., 2018; Khodasevich et
al., 2020; Menéndez-Velázquez et al., 2022). Therefore, monitoring night sky
brightness (NSB) was soon recognized as being an inevitable effort in order to
keep track of light pollution. Various organisations and individuals around
the globe started to continuously measure the NSB using different methods and
devices (Hänel et al., 2018). One of the probably most widely used devices is
the so called Sky Quality Meter (SQM). Operational SQM networks are found e.g.
in Austria (Posch et al., 2018; Puschnig et al., 2020), Spain (Zamorano et
al., 2015; Bará et al., 2019), Italy (Bertolo et al., 2019) and the
Netherlands (Schmidt & Spoelstra, 2020). Furthermore, several individuals have
mounted SQMs at various sites on Earth (e.g. Puschnig et al., 2014a, b;
Andreić, 2018). Kyba et al. (2015) have compiled many of these heterogene
datasets to study the change of light pollution on a global scale.

Early light pollution studies based on SQM data revealed the strong impact of
clouds (Kyba et al., 2011; Puschnig et al., 2014a; Jechow et al., 2019;
Ściężor, 2020) on the NSB. And more recently, the impact of several
atmospheric parameters, in particular aerosol optical depth (AOD) and
particulate matter (PM) were observed to show correlations with NSB (Ściężor &
Kubala, 2014; Posch et al., 2018; Ściężor & Czaplicka, 2020; Kocifaj &
Barentine, 2021). Also, other seasonal parameters such as e.g. surface albedo
or vegetation (Wallner & Kocifaj, 2019; Puschnig et al., 2020) were previously
discussed in the literature as potential factor causing changes in NSB.

More recently, Bara et al. (2021) and Puschnig et al. (2021) reported on the
degradation of the SQM sensitivity (i.e. darkening) with time, when the
devices are used outdoors. Using reference measurements from unused/unexposed
SQMs at the beginning and end of a multi-year time series, Bara et al. (2021)
find that readings from their SQMs, located in Galicia, need to be corrected
by approximately 0.06
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ỹr-1. Puschnig et al.
(2021) have shown that twilight may serve as a source for calibration and they
find that their readings (Vienna, Berlin and Stockholm) need to be corrected
by approximately 0.03–0.05
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ỹr-1. A combined view on
these results suggests that the “aging effect” depends on solar irradiance,
that is a function of geographic latitude.

The assessment of long-term trends of light pollution at a given site is thus
a complex problem that involves knowledge of atmospheric parameters as well as
knowledge/derivation of the instrumental darkening of the SQM (or similar
devices) with time. In this paper we aim to solve the problem via the
combination of SQM data with atmospheric parameters obtained through the
Copernicus Climate Change Service (C3S) and the Copernicus Atmosphere
Monitoring Service (CAMS; Inness et al. 2019). We establish a simple empirical
model to predict NSB variations due to changing atmospheric conditions. As a
result, the purely ALAN-driven long-term change of NSB is assessed.

The paper is organised as follows. In Section 2 we give an overview of our SQM
sites and the atmospheric parameters that are used in this study. In Section 3
we describe the data reduction routine that we apply to the SQM data to select
for moonless, clear sky measurements. The section also contains the
description of the empirical atmospheric model to predict NSB variations. The
derived long-term trends for our sites are presented in Section 4. Finally, we
discuss the results and summarize the paper in Sections 5 and 6.

##  2 Measurements and Data Products

Table 1: Basic information about our SQM measurement sites, categorized into urban, intermediate and rural ones. Code | Name | Latitude N | Longitude E | Elevation [m] | data | operating
---|---|---|---|---|---|---
|  |  |  | (above sea level) | from-to | time in yrs.
urban
STO | Stockholm (AlbaNova University Center) | N 59 21 12 | E 18 3 28 | 30 | Dec. 2014 – Dec. 2021 | 7.1
IFA | Vienna (Insitute for Astronomy) | N 48 13 54 | E 16 20 3 | 250 | Apr. 2012 – Dec. 2021 | 9.75
BA1 | Potsdam – Babelsberg | N 52 22 48 | E 13 6 22 | 90 | Jan. 2011 – Feb. 2021 | 10.2
GRA | Graz – Lustbuehel |  |  |  | Jul. 2014 – Jun. 2020 | 6.0
LSM | Linz, Schlossmuseum | N 48 18 19 | E 14 16 58 | 287 | Aug. 2014 – Dec. 2021 | 7.4
LGO | Linz, Göthestraße | N 48 18 19 | E 14 18 30 | 259 | Jan. 2014 – Dec. 2021 | 8.0
STY | Steyr | N 48 2 57 | E 14 26 32 | 307 | Aug. 2014 – Dec. 2021 | 7.4
TRA | Traun | N 48 14 8 | E 14 15 11 | 269 | Jan. 2015 – Dec. 2021 | 7.0
WEL | Wels, Rathaus | N 48 9 23 | E 14 1 29 | 317 | Aug. 2014 – Dec. 2021 | 7.4
intermediate
BRA | Braunau | N 48 15 40 | E 13 2 41 | 351 | Jan. 2016 – Dec. 2021 | 6.0
GRI | Grieskirchen | N 48 14 4 | E 13 49 33 | 336 | Jan. 2016 – Dec. 2021 | 6.0
FRE | Freistadt | N 48 30 33 | E 14 30 7 | 512 | Jan. 2016 – Dec. 2021 | 6.0
MAT | Mattighofen | N 48 5 50 | E 13 9 6 | 454 | Jan. 2016 – Dec. 2021 | 6.0
PAS | Pasching | N 48 15 31 | E 14 12 36 | 292 | Jan. 2015 – Dec. 2021 | 7.0
VOE | Vöcklabruck | N 48 0 21 | E 13 38 43 | 434 | Jan. 2016 – Dec. 2021 | 6.0
rural
FOA | Mitterschöpfl | N 15 55 24 | E 48 5 3 | 880 | Jan. 2013 – Jul. 2019 | 6.6
BOD | Nationalpark Bodinggraben | N 47 47 31 | E 14 23 38 | 641 | Jul. 2016 – Dec. 2021 | 5.5
FEU | Feuerkogel | N 47 48 57 | E 13 43 15 | 1628 | Jan. 2016 – Dec. 2021 | 6.0
GIS | Giselawarte | N 48 23 3 | E 14 15 11 | 902 | Jan. 2016 – Dec. 2019 | 4.0
GRU | Grünbach | N 48 31 50 | E 14 34 30 | 918 | Jan. 2016 – Dec. 2021 | 6.0
KID | Kirchschlag – Davidschlag | N 48 26 31 | E 14 16 26 | 813 | Jan. 2016 – Dec. 2021 | 6.0
KRI | Krippenstein | N 47 31 23 | E 13 41 36 | 2067 | Jan. 2016 – Dec. 2021 | 6.0
LOS | Losenstein, Hohe Dirn | N 47 54 22 | E 14 24 40 | 982 | Jan. 2016 – Dec. 2021 | 6.0
MUN | Münzkirchen | N 48 28 45 | E 13 33 29 | 486 | Jan. 2016 – Dec. 2021 | 6.0
ULI | Ulrichsberg, Schöneben | N 48 42 20 | E 13 56 44 | 935 | Jan. 2016 – Dec. 2021 | 6.0
ZOE | Nationalpark Zöbloden | N 47 50 18 | E 14 26 28 | 899 | Jan. 2016 – Dec. 2021 | 6.0

###  2.1 Night Sky Brightness Measurements in Stockholm, Berlin, Vienna and 23
more sites in Austria

This study is based on long-term zenithal NSB measurements obtained with Sky
Quality Meters (SQMs) located at various sites (see Table 1), including
metropolitan areas of Stockholm, Berlin and Vienna. In previous studies, we
already used parts of the data that is also included in this work, e.g.
Puschnig et al. (2014b) and Puschnig et al. (2014a) studied the influence of
the Moon, clouds and other environmental effects on the night sky brightness
over Potsdam and Vienna using the first 1–2 years of SQM data. A detailed
description and quantification of the light pollution level at the Upper
Austrian sites is found in Posch et al. (2018), who examined data obtained
during the years of 2015 and 2016.

We cover very remote locations such as Krippenstein on the Dachstein plateau
($\sim$2000 m above sea level) as well as large metropolitan areas such as
Stockholm (STO), Vienna (IFA) or Potsdam-Babelsberg (BA1), located $\sim$23 km
to the southwest of the center of Berlin.

The measurements are taken in an automated way, with LAN-attached SQM devices
(model SQM-LE) located in weather-proof housings. The SQMs situated in Upper
Austria are run by the provincial government of Upper Austria. They take NSB
measurements every minute. Most other SQMs (STO, IFA, FOA, GRA) provide a
reading every 7 seconds which corresponds to a frequency of 0.143 Hz and BA1
even takes a measurement every $\sim$2 seconds.

###  2.2 Atmospheric data products

We make use of open access climate variables from ERA5 (Hersbach et al.,
2018), the fifth major global reanalysis data produced by the European Centre
for Medium-Range Weather Forecasts (ECMWF). ERA5 is developed through the
Copernicus Climate Change Service (C3S). The data are based on a reanalysis of
a large set of ground-, air- and satellite-based measurements. Data are
available with hourly validity time at a spatial resolution of 0.28° x 0.28°
in latitude and longitude, corresponding to $\sim$30 km x 30 km).

Using the python package cdsapi provided by the Atmosphere Data
Store111https://ads.atmosphere.copernicus.eu, we downloaded for all our sites
and dates more than 70 available atmospheric parameters. From those we
initially identified quantities that may have impact on zenithal NSB
measurements. This pre-selection is mainly driven by our understanding of how
the environment and the atmosphere impact the NSB, i.e. we chose parameters
that may alter the NSB via scattering of light (e.g. aerosols and particles)
as well as parameters that may enhance the fraction of upward light (e.g.
albedo, snow cover) or enhance the fraction of light that is reflected back to
the ground (e.g. cloud cover). An overview of the selected parameters is found
in Table 2.

We note that aerosol optical depths and particulate matter are not available
from ERA, but are provided through CAMS (Inness et al., 2019), the Copernicus
Atmosphere Monitoring Service global reanalysis. Daily Data products are
available in steps of six hours. Due to a bug that caused AOD and PM to not be
computed for analysis
time222https://confluence.ecmwf.int/pages/viewpage.action?pageId=153393473, we
had to use the 3-hour forecasts with a validity time of 3am UTC. CAMS native
spatial resolution is 0.8° before 21 June 2016, and 0.4° henceforth,
corresponding to $\sim$80 km x 80 km and $\sim$40 km x 40 km.

####  2.2.1 Description of selected atmospheric parameters

An overview of our initial set of parameters used for statistical analysis is
given in Table 2. In the following we briefly describe the physical
underpinning of the most relevant parameters.

Albedo (aluvd) is defined as the fraction of incident radiation that is
reflected by a surface. Its numerical value thus ranges from 0 (no reflection)
to 1 (all incident radiation is reflected). It varies with the type of surface
(e.g. its roughness) and wavelength. For example, the broad-band albedo of
grassland is a few percent only (Briegleb & Ramanathan, 1982; Briegleb et al.,
1986; Coakley, 2003), with twice as much reflection in the near infrared than
in the visible spectral range. On the other side, surfaces covered with snow
may reflect more then 90 percent of the incident radiation, with snow albedo
being significantly higher in the visible range than in the near infrared
(Roesch et al., 2002). In this study, we make use of the UV-optical albedo for
diffuse radiation (aluvd), that is measured within a range of 300 to 700nm and
thus matches the SQM band very well. It is thus expected that variations in
surface albedo positively correlate with the night sky brightness.

Leaf Area Index (lai) is defined as the one-sided green leaf area per unit
ground area (Boussetta et al., 2011). It is thus a dimensionless number
($\mathsf{m^{2}\ m^{-2}}$). In ERA5, the leaf area index is split into high
(laihv) and low vegetation (lailv). High vegetation consists of evergreen
trees, deciduous trees, mixed forest/woodland, and interrupted forest and low
vegetation covers grass, shrubs as well as water and land mixtures. For our
study we use the sum of both available quantities. High vegetation impacts the
NSB via blocking of upward light, while low vegetation has an impact on the
reflection of downward light. We expect a negative correlation between NSB and
the leaf area index.

Both, the leaf area index and albedo are based on observations with the
VEGETATION sensor on board the SPOT satellite, that carries out measurements
every 10 days using a composite observation from a 30 days moving window at
1/112 degree spatial resolution, corresponding to roughly 1 km at the equator
(Boussetta et al., 2014).

Black Carbon (bcaod550) is a direct consequence of anthropogenic activities,
in particular black carbon is a result of incomplete combustion of fossil
fuels (38%), biomass (42%) and biofuels (20%) (Bond et al., 2004; Cao et al.,
2014; Schwarz et al., 2006). Black carbon in the atmosphere is capable of
absorbing and scattering of radiation. The atmospheric black carbon aerosol
optical depth measured at 550 nm (bcaod550) has thus the potential of having a
large impact on our NSB measurements. Note that black carbon may also have an
impact on the surface albedo, e.g. when it is deposited on snow and ice.

Table 2: Overview of atmospheric parameters, i.e the features in the model aluvd | UV-optical albedo for diffuse radiation
---|---
lai | leaf area index
bcaod550 | black carbon aerosol optical depth
duaod550 | dust aerosol optical depth
ssaod550 | sea salt aerosol optical depth
omaod550 | organic matter aerosol optical depth
pm10 | particulate matter 10$\mu$m
pm2p5 | particulate matter 2.5$\mu$m
pm1 | particulate matter 1$\mu$m
tcwv | total column water vapour
tcw | total column water
wind | wind 10m
tco3 | total column ozone
tcc | total cloud cover
sd | snow depth

##  3 Methods

###  3.1 Atmospheric model

We aim to find a set of atmospheric parameters that most directly impact
(zenithal) NSB measurements. In particular, we attempt to separate fundamental
correlations between the NSB and the atmosphere from those that are indirect
consequences of covariance among atmospheric metrics. We distinguish these
underlying relations through variable selection: For each NSB data (as a
target variable), we compose an empirical predictive model using a set of
atmospheric parameters (feature variables) that carry most predictive power.
This is an effective way to collapse a high-dimensional data set into a
concise model.

The basis of this analysis is a multi-variable penalized linear regression.
That is, we restrict the model functional forms to simple linear combinations
of variables (including an intercept term). The regression is done
independently for each NSB measurement/chunk ($\Delta_{\mathrm{sky-obs}}$) as
a target variable, using all the atmospheric parameters plus time (to account
for temporal trends) as available features. With this regression setup, we
perform a lasso model fit (Tibshirani, 1996) and use the Bayesian Information
Criterion (BIC; Schwarz 1978) for automated feature/model selection. This is
implemented with the LassoLarsIC function in the scikit-learn Python package.
In detail, for a linear predictive model with the form
$\hat{y}_{i}=\beta_{0}+\sum_{j=1}^{m}\beta_{j}x_{ij}$, the lasso regression
minimizes the following function:

$\frac{1}{2n}\sum_{i=1}^{n}\ (y_{i}-\hat{y}_{i})^{2}\ +\alpha\ \sum_{j=1}^{m}\
|\beta_{j}|.$ (1)

The indices i and j run through the NSB data and the atmospheric features
respectively. The parameter $\alpha$ is a hyper-parameter, so that the second
term in the equation adds a penalty for the use of any non-zero slope in the
fitted model. This particular “regularization” term is the reason that the
lasso as a regression method can also be used for variable selection. The
lasso regression yields a best-fit model that minimizes the equation for each
choice of the $\alpha$ parameter.

In practice, we produce 16 different predictive atmospheric models. Starting
with a model based on all 15 atmospheric and ground features plus time, we
expect to find the best agreement between observed NSB variations and the
models. However, we also aim to quantify the impact of each parameter. Hence,
we iteratively remove input features and re-calculate the models. This will
allow us to use the residual after subtracting the model from the data as a
measure of the importance of individual input features for the NSB modeling.

Note that due to the “regularization” term not all input features may survive
the lasso method and only those parameters that have most direct impact on the
NSB will be included in the final model.

Figure 1: Measuring the SQM aging effect from twilight observations. Left:
Germany, Potsdam-Babelsberg (BA1; 10.2 years) Right: Austria, Bodinggraben
(BOD; 5.5 years). Topmost panels: Patat model minus observed clear-sky SQM
twilight data on y-axis versus time on x-axis. A strong degradation with time
is revealed. Each point is a 5-minute-mean (BA1) or 10-minute-mean (BOD)
brightness. 2nd panels from top: Prediction from empirical atmospheric model
plus linear time trend. 3rd panels from top: The residual after subtraction of
the top-panel data from the model in the 2nd panel. The red line and area show
a polynomial fit and its standard deviation respectively. Bottom panels: Final
scatter after subtraction of the atmospheric model and the linear and
polynomial trends.

###  3.2 Correcting SQM readings for degrading sensitivity with time using the
twilight method

We apply an updated version of the “twilight” calibration as described in
Puschnig et al. (2021) in order to correct our multiple-year lasting SQM
measurements for potential temporal changes in detector sensitivity. This is
the “aging effect”, which is probably caused by changing transmission of the
SQM housing window due to UV light exposure over several years. In brief, we
utilize SQM observations of the zenithal NSB and compare the measurements to
the twilight model of Patat et al. (2006). Moon and Sun altitudes are
calculated using the Python ephem package, which provides an accuracy of
approximately 1 arcsec. Note that one could also use an inter-comparison of
SQM data obtained during twilight when the sun was at the same altitude.
Assuming that the twilight sky brightness is not affected by ALAN and remains
constant, one can reveal any underlying temporal sensitivity change of the
measurement system. In addition to this procedure, we now also account for
atmospheric changes that may impact the twilight observations. Thus, we
utilize the same atmospheric model as described in Section 3.1. Any linear
temporal trend obtained from the modeling procedure is then attributed to the
“aging effect”.

Since we recognized at some stations (e.g. IFA) remaining non-linear trends,
we additionally fit polynomials up to 3rd order to the residuals. If the
subtraction of a polynomial further reduces the scatter (by more than 0.5
percent), we proceed and use also the polynomial as a correction function,
accounting for non-linear effects.

A quantification of the degradation of SQM sensitivity with time is seen in
Figure 1 for Potsdam-Babelsberg near Berlin and a rural site in Austria
(Bodinggraben). Note the strong seasonal variation of NSB during twilight and
how accurate the atmospheric model is able to predict them for both the urban
site and the rural site. The histogram in Figure 2 shows the distribution of
linear aging slopes found for all our 26 SQM stations. We find a relatively
large variation between individual SQM sites with aging slopes ranging from
zero to $\sim$-0.075 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$
yr-1. On average, the aging effect leads to a darkening of
$\sim$-0.031$\pm$0.020 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$
yr-1.

Figure 2: Histogram of the linear part of the aging slopes derived using the
“twilight” method for our SQM sites. The vertical black line denotes the mean
slope. Note that the unit on the x-axis is
m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$.

###  3.3 Extraction of moonless, clear sky SQM measurements

We are interested in the anthropogenic contribution to the NSB. Thus, we
compare our SQM measurements to a modeled zenithal NSB calculated from sky
spectra available through the SKYCALC sky model (Noll et al., 2012; Jones et
al., 2013). The whole procedure of how we extracted SQM magnitudes from the
model spectra is described in the Appendix A. From this point on, we proceed
with the difference between the sky model (sky) and the observed NSB (obs):
($\Delta_{\mathrm{sky-obs}}$). Note that – given our model constraints – the
SKYCALC sky models lead to yearly peak and valley zenithal NSBs (due to
changing starlight+zodiacal light) of 21.63 and 21.87
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ for our sites. Comparing
our yearly peak-to-valley difference of 0.24
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ to the “GAIA map of the
brightness of the natural sky” (Masana et al., 2022), shows that this is a
relatively low value. Figure 3 of Masana et al. (2022) suggests a maximum
yearly V-band variations of approx. 0.6 mag arcsec-2. This discrepancy arises
due to several facts. First, we evalute the SKYCALC models using (for the sake
of simplicity) constraints such as a fixed precipitable water vapor value of 5
mm rather than a time-dependet value. This has impact on the scattering of
light in the atmosphere and will lead to lower NSB variations over the year.
Second, our sites are located at geographic latitudes between 50 and 60
degrees. This is higher than the example shown in Figure 3 of Masana et al.
Thus the contribution of zodiacal light to the zenithal NSB becomes lower at
our sites (Masana et al., 2020, see Figure 10), further reducing the overall
yearly zenithal NSB variation.

Our data reduction routine starts with splitting the SQM data into chunks of
45 minute length. We have chosen this time span to be short enough to allow
for multiple sampling points each night (even during summer, except for
Stockholm) and to be long/large enough to avoid stochasticity. For each of
these time chunks we calculate the mean NSB and the standard deviation. Also a
linear fit is performed and the maximum deviation from the fit line is
determined. For the mean time of each chunk we further lookup the moon
altitude. That way, we are prepared to downselect for observations obtained
during moonless and clear sky conditions. We only keep chunks that fulfill the
following constraints: i) the Moon is below the horizon, ii) the maximum
deviation from the linear fit line of any single measurement is lower than
0.04 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$, iii) the standard
deviation of the measurements within the chunk is less than 0.02
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ for SQMs with a high
sampling rate (BA1, IFA, STO, FOA) and 0.06
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ for SQMs with a sampling
frequency of only 1/minute, iv) the slope of the linear fit is shallower than
0.13 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ per hour.

This procedure basically reduces the SQM data to chunks of 45 minute length
within which the NSB remains almost constant. This is the case only when no
clouds are present in zenith. Note that the capability of using the SQM as
cloud detector was previously demonstrated by Cavazzani et al. (2020).
Finally, we lookup meteorological parameters for each of the data chunks and
reject data chunks where i) the large-scale total cloud cover was larger than
50 percent and ii) the snow depth was larger than 5 cm (avoiding SQM readings
when the sensor was potentially covered with snow).

Figure 3: Left: Germany, Potsdam-Babelsberg (BA1). Right: Austria,
Bodinggraben (BOD). Top panels: The difference between the predicted (skycalc
model) and the observed clear-sky night sky brightness data within the SQM
sensitivity range as a function of time. Each point represents the mean NSB
within time chunks of 45 minute length. The gray area shows the
$\pm$1.5$\upsigma$ range of a linear least square fit to the points. Mid
panels: The best-fit (penalized linear regression) model based on 15
atmospheric parameters obtained via the Copernicus Climate Change Service
(C3S) Data Store. Bottom panels: The residual after subtraction (model minus
data). The panels to the right of the time series are histograms of the data
points, using a binwidth of
0.05${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$. Note the observed
seasonal variations of NSB for the rural station in Austria and how accurate
the atmospheric model is able to predict them.

###  3.4 Error estimation

In order to estimate the error of our final trend slopes expressed in
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1 we take into
account uncertainties due to the aging correction and the trend fitting. Both
errors are estimated through examination of linear fits through the residuals,
e.g. bottom panels in Figures 1 and 3. We calculate the linear regression
uncertainties from the mean squared errors which are the trace (equal to the
sum of elements on main diagonal) of the error covariance matrix. That way, we
find 1-sigma uncertainties due to the aging correction ranging from 0.001 to
0.006 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1 with a mean
uncertainty of 0.0035 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$
yr-1. The 1-sigma uncertainties resulting from the linear trend analysis vary
between 0.001 and 0.008 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$
yr-1 with a mean uncertainty of 0.0045
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1. Using a
conservative approach, we finally calculate the absolute maximum 1-sigma error
for our routine simply via addition: 0.0035 + 0.0045 = 0.0080
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1. We thus conclude
that our method is capable of detecting trends that are shallower/steeper than
$\pm$0.016 ${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1,
corresponding to the $\pm$2-sigma level. Expressed on a linear scale, our
trend detection limit is thus $\pm$1.5 percent per year.

##  4 Results

###  4.1 Long-term trends

After selecting aging-corrected SQM data chunks as explained in Section 3, one
can easily see in the top panel of Figure 3 that even under clear and moonless
conditions individual SQM measurements are prone to large variations with
peak-to-valley differences up to
$\sim$1${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$, corresponding to
a factor of 2.5. At rural sites such as e.g. BOD, the scatter is dominated by
a seasonal variation, caused by variations of albedo (e.g. enhancing NSB
during winters due to snow cover) and vegetation (darkening during summers due
to blocking of light on leafs). In urban areas (e.g. BA1) such seasonal
effects are typically weaker and other parameters seem to gain importance,
e.g. atmospheric aerosols or particulate matter.

Pronounced variations and the scattering of data points hamper the assessment
of long-term trends. For example, non-uniform, stochastic sampling of a
periodic function (seasonal variations) may introduce a bias that leads to a
spurious trend when performing a linear regression. For that reason, we aim to
model and remove the impact of the atmosphere. To do so, we use atmospheric
products freely available from the Copernicus Earth observation program. Our
empirical atmospheric model is then found from a penalized linear regression
method (lasso) as explained in Section 3.

The middle panel in Figure 3 makes evident that our model is capable of
predicting the bulk of the seasonal NSB variations solely from the atmospheric
data. The long-term linear trends we find for Potsdam-Babelsberg (BA1) and
Bodinggraben (BOD) are 45$\pm$16 and 27$\pm$16
m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1, respectively.
This corresponds to an increase of $\sim$4$\pm$1.5% and 2$\pm$1.5% per year.

The average increase in light pollution at our 11 rural sites is 1.7$\pm$1.5%
per year. At our nine urban sites we measure an increase of 1.8$\pm$1.5% per
year and for the remaining six intermediate sites we find an average increase
of 3.7$\pm$1.5% per year. These numbers correspond to doubling times of 41, 39
and 19 years. Results for all 26 stations are found in Table 3.

Table 3: For each SQM station in column 1, the linear slope of the aging effect given in m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ year-1 is shown in column 2. Column 3 contains the order of the polynomial used to account for non-linear sensitivity changes of the SQM. The long-term linear trend slope in m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ year-1 is found in column 4. All features that survived the lasso regression are listed in column 5 and the total number of data chunks are given in column 6. The last column contains the Pearson R correlation coefficient for the full model (using 16 features). The best fits with R>0.8 are BOD, FEU, KRI, STO. Code | aging | polyn. | trend | features that | data | R
---|---|---|---|---|---|---
| slope | order | slope | survived lasso | count |
urban
BA1 | -48 | 2 | 45 | time, lai, bcaod550, duaod550, omaod550, ssaod550 | 1819 | 0.74
GRA | -73 | 3 | 49 | time, aluvd, lai, bcaod550, pm1, duaod550, omaod550, ssaod550, pm2p5, tcw, sd, wind10, tco3, tcc | 846 | 0.58
IFA | -68 | 2 | 35 | time, aluvd, lai, bcaod550, pm1, duaod550, omaod550, pm10, pm2p5, tcwv, wind10, tco3, tcc | 1506 | 0.59
LGO | 6 | 3 | -12 | aluvd, lai, pm1, duaod550 | 315 | 0.59
LSM | -35 | 3 | -16 | aluvd, lai, sd | 424 | 0.6
STY | -15 | 2 | 19 | time, aluvd, lai, bcaod550, duaod550, omaod550, ssaod550, pm2p5, sd, wind10, tco3 | 1335 | 0.51
TRA | -37 | 3 | 29 | time, aluvd, lai, duaod550, omaod550, ssaod550, wind10, tco3 | 410 | 0.53
WEL | -29 | 3 | -9 | time, aluvd, lai, pm1, duaod550, omaod550, ssaod550, wind10 | 1101 | 0.64
STO | -43 | 0 | 37 | time, aluvd, lai, duaod550, omaod550, tcwv, sd | 913 | 0.82
intermediate
BRA | -45 | 3 | 31 | time, aluvd, lai, pm1, duaod550, omaod550, ssaod550, tco3 | 970 | 0.58
FRE | -36 | 0 | 21 | time, aluvd, lai, bcaod550, duaod550, omaod550, ssaod550, sd, tco3 | 960 | 0.61
GRI | -45 | 0 | 46 | time, aluvd, lai, bcaod550, pm1, ssaod550, pm10, pm2p5, tcwv, sd, wind10, tco3 | 893 | 0.67
MAT | -48 | 3 | 90 | time, aluvd, lai, bcaod550, duaod550, ssaod550, tcwv, sd, wind10, tco3, tcc | 979 | 0.64
PAS | -29 | 3 | 7 | time, aluvd, lai, bcaod550, pm1, duaod550, omaod550, pm10, pm2p5, tcwv, sd, wind10, tco3, tcc | 564 | 0.6
VOE | -50 | 0 | 50 | time, aluvd, lai, bcaod550, duaod550, omaod550, ssaod550, pm10, tcwv, sd, wind10, tco3, tcc | 678 | 0.54
rural
FOA | -28 | 3 | 7 | time, aluvd, lai, bcaod550, duaod550, omaod550, ssaod550, pm10, tcw, sd, wind10, tco3, tcc | 1122 | 0.65
FEU | -5 | 3 | 21 | aluvd, tco3 | 848 | 0.81
GIS | 0 | 2 | 62 | time, aluvd, lai, bcaod550, duaod550, pm10, wind10, tco3, tcc | 343 | 0.59
GRU | -17 | 0 | 11 | aluvd, lai, bcaod550, duaod550, omaod550, tcwv, sd, wind10, tco3, tcc | 1265 | 0.63
KID | -11 | 0 | -11 | time, aluvd, lai, bcaod550, duaod550, omaod550, tcwv, sd, tco3, tcc | 1069 | 0.55
KRI | -24 | 3 | 22 | time, aluvd, lai, bcaod550, tco3 | 800 | 0.82
LOS | -51 | 0 | 1 | aluvd, lai, omaod550, ssaod550, tcwv, tco3 | 913 | 0.54
MUN | -58 | 3 | 61 | time, aluvd, lai, bcaod550, duaod550, sd, tco3 | 860 | 0.72
BOD | -31 | 3 | 27 | time, aluvd, lai, bcaod550, omaod550, ssaod550, tco3 | 874 | 0.84
ZOE | -17 | 3 | 5 | aluvd, lai, duaod550, omaod550, tco3 | 817 | 0.78
ULI | -1 | 3 | 3 | time, aluvd, lai, bcaod550, duaod550, omaod550, pm10, tcw, sd, wind10, tco3, tcc | 908 | 0.68

###  4.2 Impact of the atmosphere

As explained in Section 3, we perform a multi-variate penalized linear
regression to identify the importance of individual atmospheric parameters. In
detail, we calculate models with increasing complexity, i.e. with an
increasing number of atmospheric input features. Starting with a model based
on only a single variable (time), any potential long-term linear trend is
revealed. Then, we consecutively add single atmospheric features (those listed
in Table 2) to the model in order to study the relative importance of each
parameter. That is, we compare the scatter (standard deviation) of the
residuals after subtracting the model from the data.

As shown in Figure 4, additional features reduce the residual scatter. In
particular, the Figure shows that at all our sites, the leaf area index (lai)
and surface albedo (aluvd) are the most important model input parameters that
reduce the final scatter by a large fraction. Snow depth (sd) in principle may
firmly enhance the NSB. However, we did not expect that sd has a strong impact
on our models, since for the bulk of our locations sd is typically zero most
of the time. Moreover, we rejected data points when the snow cover was more
than 5 cm (to avoid observations when the SQM sensor was covered with snow).
Yet, we do see a decreasing residual scatter due to sd for some rural sites
(BOD, FEU, KRI, ZOE) and we observe that snow depth even has a strong effect
in Stockholm (STO). Variations of albedo mostly impact our rural sites, while
urban and intermediate sites are less affected. In contrast to that, the leaf
area index (lai) helps in reducing the scatter of all sites by a large degree,
be it rural, intermediate or urban.

Atmospheric black carbon (bcaod550) is an important ingredient for our models
of urban and (in particular) intermediate sites, while it has almost no impact
on rural sites. Interestingly, our coarsly spatially resolved atmospheric
data, does not suggest any strong correlation between NSB and large
particulate matter (pm10). However, for smaller grains, we do see that pm2p5
and pm1 reduce the residual scatter of the bulk of our urban and intermediate
sites. Organic matter in the atmosphere (omaod550) plays mostly only a minor
role at all rural sites, but at one station (GIS), omaod550 seems to be
important for the modeling.

Ozone makes no difference for the modeling of our urban and intermediate
sites. However, models that include ozone as a feature reduce the residual
scatter of some rural sites such as ULI, GIS, GRU. Finally, we find that the
10-meter wind speed (wind10) improves the match between the models and NSB
observations for some of our sites (GRA, GRI, GIS). The correlation between
wind speed and NSB is a non-causal secondary effect, that may be the result of
an increase in aerosol transport when wind speed increases. We reckon that due
to the coarser temporal resolution of our AOD data (one measurement per six
hours) compared to the wind data (hourly data), at sites where aerosol
abundance correlates with wind speed, the wind parameter may serve as a proxy
for AOD variations at shorter time scales.

For each of our models and sites we calculate the Pearson correlation
coefficient $R$ (see Table 3). While the correlation between the SQM data and
a solely linear trend is weak ($R\sim 0.3$ on average), the correlation
between the full atmospheric model and the SQM observations becomes strong
($R\sim 0.7$ on average). For some of the sites (BOD, FEU, KRI, STO) the
Pearson correlation coefficient even increases to values greater than 0.8,
indicative for a very strong correlation.

Figure 4: Relative impact of atmospheric parameters. The number of model input
features increases from left to right along the x-axis, starting with a solely
linear time dependence (very left; labelled as no model) to 16 features for
the full model (right end; time plus 15 atmospheric features). The y-axis is
the residual scatter (standard deviation) found after subtracting the model
from the data. The thick curves with points indicate the mean scatter over all
stations in the panel. The numbers in the boxes are the Pearson correlation
coefficients (R) without any atmospheric model (left) and with the full model
(right).

##  5 Discussion

Our long-term analysis of zenithal NSB measurements reveal that it is not the
urban regions that show the largest increase in light pollution, but rather
the intermediate regions. This may be explained by the fact that the
installation of additional lighting points in urban areas contributes
relatively less to the NSB, as in urban areas the overall brightness level is
already extremely high. On the other side, when intermediate areas develop,
new or upgraded lights have a much larger impact on the NSB, as it is
relatively darker than in urban areas. More so, in rural regions. There, even
only few additional lights may have a recognizable impact. Having said that,
Stockholm is probably an exception of that average rule, because it is our
brightest site and still shows a very strong increase.

For three of our stations (BA1, BOD, MAT) we utilize the “Radiance Light
Trends” web application333https://lighttrends.lightpollutionmap.info, which
provides monthly mean radiance measurements obtained with the “Day-Night Band”
(DNB) of the “Visible Infrared Imaging Radiometer Suite” (VIIRS) aboard the
Suomi-NPP satellite. The resulting long-term radiance trends are shown in
Figure 5. The VIIRS monthly averages show a large scatter and no seasonal
variation. This is caused by several factors. First of all, VIIRS data is only
available for our sites between September and March, because only during these
months the satellite passes during astronomical night. This may drastically
reduce the dynamic range of the observations and basically filters out
seasonal trends. Thus, our VIIRS data contain only measurements of the winter
half-year, which are more prone to large changes in albedo. As revealed by our
analysis, albedo has a strong impact on the NSB. Thus, data obtained during
the winter-half-year tend to show larger scatter. Moreover, the variation of
the satellite’s viewing angle adds another factor that increases the scatter.
The monthly averaging may reduce some of these problems (e.g. the viewing
angle), but still large uncertainties as seen in Figure 5 remain, making a
long-term trend analysis challenging. On top of that, the VIIRS/DNB spectral
response is very different from the SQM. While the SQM is most sensitive
between 400 and 600nm, the VIIRS/DNB passband covers a range between 500 and
900nm, which is most problematic in case of a change of the emitted spectral
energy distribution towards blue colors (e.g. upgrade to 4000K LEDs). This
should be kept in mind when comparing our SQM-based results to those from
VIIRS. Therefore, we do not expect to see a perfect match between trends
obtained from such different instruments. However, we do expect that direction
and magnitude of the trends are comparable. After conversion to magnitudes,
the VIIRS-based trend slopes are 29$\pm$30, 133$\pm$143 and 94$\pm$21 mmag
arcsec-2 yr-1 for BA1, BOD and MAT respectively, i.e. the uncertainties are
quite large in relation to our SQM-based study. Only for one of the stations,
MAT, a comparison can be made. Interestingly, in this case, the results do
match very well. From VIIRS an increase of 94$\pm$21 mmag arcsec-2 yr-1 (over
the DNB bandpass) is measured, while our SQM analysis suggests a slope of
90$\pm$16 m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1. This
may be an indication that the color of the night sky at this site has not
changed much.

Finally, we briefly discuss the SQM aging or time-dependent loss of
sensitivity as deduced from the twilight method. Figure 2 shows that the
correction slopes range from zero to almost $-$80
m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1 with a mean and
standard deviation of $-$31$\pm$20
m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1. It was previously
speculated by Bara et al. (2021) and Puschnig et al. (2021) that the darkening
may be related to the total (time-dependent) amount of (UV) light exposure at
a given site. The scatter of the aging slopes among our stations, however,
does not support such a simple picture. Some other factors, probably dust and
polls, may play a crucial role. Regardless of what causes the loss of
sensitivity, we are confident that the post-calibration procedure via twilight
measurements is capable of correcting the change in sensitivity. An example
that provides evidence for the robustness of our method is the comparison of
the results obtained at the stations LGO and LSM. These two SQMs are situated
in the same city (Linz), and are thus only $\sim$1km apart. While for one of
the sites (LGO) the aging effect is even below the recognition limit, for the
other one (LSM), an aging slope of $-$35
m${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ yr-1 is found. While
the correction slopes are quite different, the final trends are in good
agreement. A comparison of several SQMs operating in parallel at different
sites within a small area could potentially provide more insight into the
aging effect in the future.

Figure 5: Radiance trends for three of our sites obtained through the
“Radiance Light Trends” web application. A radius of 10km around our SQM sites
was chosen for summing up radiances. Each red point is a monthly average and
the gray shaded area represents the uncertainty of a linear fit.

##  6 Summary and Conclusion

We analysed long-term (4-10 years) SQM measurements obtained at 26 different
sites, of which 24 are located in Austria, one in Stockholm (Sweden) and one
in Potsdam-Babelsberg (Germany). We utilize the “twilight method” to correct
for changing detector sensitivity with time, i.e. the “aging effect”. Using
aging-corrected data obtained under clear and moonless sky, long-term trends
of light pollution are derived for our sites. Categorizing our sites into
rural, urban and intermediate ones, we find an average increase in light
pollution of approximately 1.7, 1,8 and 3.7 percent per year respectively,
with an estimated $\pm$2-sigma uncertainty of $\pm$1.5% per year. The
corresponding doubling times are 41, 39 and 19 years.

Furthermore, we establish an empirical atmospheric model which allows us to
investigate the relative impact of 15 different atmospheric parameters on the
night sky brightness. We find that surface albedo and vegetation have by far
the largest impact on the zenithal night sky brightness. Additionally, black
carbon and organic matter aerosols are important at urban and rural sites,
respectively. Snow depth was found to be important for some sites, while the
total column of ozone leaves impact on some rural places.

In the paper we have shown that large-scale (30–40 km resolution) atmospheric
parameters obtained through the Copernicus Climate Change Service and the
Copernicus Atmosphere Monitoring Service may serve as a basis of a predictive
model to constrain the (zenithal) NSB. However, we also recognized that our
empirical model correlates better with SQM observations obtained at rural
sites, i.e. where albedo and vegetation have strongest impact, while aerosols
and particulate matter play only a minor role. Obviously, at urban sites the
situation is different. There, aerosols and particulate matter are typically
more abundant, and vary on shorter times and scales. The weaker correlation
between our model and the observations at urban sites is likely due to the
relatively low spatial resolution of the available parameters (aerosol optical
depths and particulate matter). Matching of high-quality aerosol optical
depths and particulate matter would most likely lead to improved predictions
also at our urban sites. This task remains to be done in the future.

## Acknowledgements

We thank the referee for carefully reading our manuscript and for constructive
comments and suggestion. We acknowledge support from the Alva and Lennart
Dahlmark research grants of the Department of Astronomy at Stockholm
University, through which one of the SQMs (STO) and its housing were financed.
We are further thankful to Sergio Gelato (Univ. of Stockholm), who helped with
the setup of the SQM in Stockholm (hard- and software).

Part of this work has been generated using Copernicus Climate Change Service
information (2022) This research made use of SciPy (Jones et al., 2001) and
NumPy (Van Der Walt et al., 2011).

## Data Availability

The SQM data underlying this article are available for download through
repositories at https://www.land-oberoesterreich.gv.at/115999.htm and
https://astro.univie.ac.at/en/science-communication/reading-material/light-
pollution/. Our reduced SQM data can be made available upon reasonable
request.

The third party meteorological data used in this article are available through
the Climate Data Store Application Program Interface (CDS API) using Python
(e.g. pip install cdsapi).

## References

  * Andreić (2018) Andreić v., 2018, Night sky brightness above Zagreb 2012.-2017, doi:10.48550/ARXIV.1803.00299, https://arxiv.org/abs/1803.00299
  * Bará et al. (2019) Bará S., Lima R. C., Zamorano J., 2019, Sustainability, 11
  * Bara et al. (2021) Bara S., Marco E., Ríbas S., Garcia M., Sanchez de Miguel A., Zamorano J., 2021, International Journal of Sustainable Lighting, 23, 1
  * Bertolo et al. (2019) Bertolo A., Binotto R., Ortolani S., Sapienza S., 2019, Journal of Imaging, 5
  * Bessell (1990) Bessell M. S., 1990, Publications of the Astronomical Society of the Pacific, 102, 1181
  * Bessell (2005) Bessell M. S., 2005, Annual Review of Astronomy and Astrophysics, 43, 293
  * Bond et al. (2004) Bond T. C., Streets D. G., Yarber K. F., Nelson S. M., Woo J.-H., Klimont Z., 2004, Journal of Geophysical Research: Atmospheres, 109
  * Boussetta et al. (2011) Boussetta S., Balsamo G., Beljaars A., Kral T., Jarlan L., 2011, Impact of a satellite-derived Leaf Area Index monthly climatology in a global Numerical Weather Prediction model, doi:10.21957/h7n0ilfkp, https://www.ecmwf.int/node/8339
  * Boussetta et al. (2014) Boussetta S., Balsamo G., Dutra E., Beljaars A., Albergel C., 2014, Analysis of surface albedo and Leaf Area Index from satellite observations and their impact on numerical weather prediction, doi:10.21957/otwcakuu3, https://www.ecmwf.int/node/12032
  * Briegleb & Ramanathan (1982) Briegleb B., Ramanathan V., 1982, Journal of Applied Meteorology, 21, 1160
  * Briegleb et al. (1986) Briegleb B. P., Minnis P., Ramanathan V., Harrison E., 1986, Journal of Climate and Applied Meteorology, 25, 214
  * Cao et al. (2014) Cao X., Liang J., Tian P., Zhang L., Quan X., Liu W., 2014, Atmospheric Pollution Research, 5, 601
  * Cavazzani et al. (2020) Cavazzani S., Ortolani S., Bertolo A., Binotto R., Fiorentin P., Carraro G., Saviane I., Zitelli V., 2020, Monthly Notices of the Royal Astronomical Society, 493, 2463
  * Chepesiuk (2009) Chepesiuk R., 2009, Environmental health perspectives, 117, A20
  * Cho et al. (2015) Cho Y., Ryu S.-H., Lee B., Kim K., Lee E., Choi J., 2015, Chronobiology international, 32, 1
  * Cinzano (2005) Cinzano P., 2005, Technical report, Night Sky Photometry with Sky Quality Meter, http://www.inquinamentoluminoso.it/download/sqmreport.pdf. Dipartimento di Astronomia, Vicolo dell Osservatorio 2, I-35100 Padova, Italy, Istituto di Scienza e Tecnologia dell Inquinamento Luminoso, Via Roma 13, I-36106 Thiene, Italy, http://www.inquinamentoluminoso.it/download/sqmreport.pdf
  * Coakley (2003) Coakley J., 2003, in Holton J. R., ed., , Encyclopedia of Atmospheric Sciences. Academic Press, Oxford, pp 1914 – 1923, doi:https://doi.org/10.1016/B0-12-227090-8/00069-5, http://www.sciencedirect.com/science/article/pii/B0122270908000695
  * Eisenbeis (2006) Eisenbeis G., 2006, Ecological consequences of artificial night lighting. Island Press, pp 281–304
  * Garcia-Saenz et al. (2018) Garcia-Saenz A., et al., 2018, Environmental Health Perspectives, 126, 047011
  * Haim & Portnov (2013) Haim A., Portnov B. A., 2013, Light Pollution as a New Risk Factor for Human Breast and Prostate Cancers, 1st ed. 2013. edn. Springer Netherlands
  * Hänel et al. (2018) Hänel A., et al., 2018, Journal of Quantitative Spectroscopy and Radiative Transfer, 205, 278
  * Hersbach et al. (2018) Hersbach H., et al., 2018, Operational global reanalysis: progress, future directions and synergies with NWP, doi:10.21957/tkic6g3wm, https://www.ecmwf.int/node/18765
  * Hölker et al. (2010) Hölker F., Wolter C., Perkin E., 2010, Trends in ecology & evolution, 25, 681
  * Inness et al. (2019) Inness A., et al., 2019, Atmos. Chem. Phys., 19, 3515
  * Jechow et al. (2019) Jechow A., Hölker F., Kyba C. C. M., 2019, Scientific Reports, 9, 1391
  * Jones et al. (2001) Jones E., Oliphant T., Peterson P., Others 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Jones et al. (2013) Jones A., Noll S., Kausch W., Szyszka C., Kimeswenger S., 2013, A&A, 560, A91
  * Khodasevich et al. (2020) Khodasevich D., Tsui S., Keung D., Skene D. J., Martinez M. E., 2020, medRxiv
  * Kocifaj & Barentine (2021) Kocifaj M., Barentine J. C., 2021, Scientific Reports, 11, 14622
  * Kyba et al. (2011) Kyba C. C. M., Ruhtz T., Fischer J., Hölker F., 2011, PLOS ONE, 6, 1
  * Kyba et al. (2015) Kyba C. C. M., et al., 2015, Scientific Reports, 5
  * Longcore & Rich (2004) Longcore T., Rich C., 2004, Frontiers in Ecology and the Environment, 2, 191
  * Masana et al. (2020) Masana E., Carrasco J. M., Bará S., Ribas S. J., 2020, Monthly Notices of the Royal Astronomical Society, 501, 5443
  * Masana et al. (2022) Masana E., Bará S., Carrasco J. M., Ribas S. J., 2022, International Journal of Sustainable Lighting, 24
  * Mathews et al. (2015) Mathews F., Roche N., Aughney T., Jones N., Day J., Baker J., Langton S., 2015, Philosophical Transactions of the Royal Society B: Biological Sciences, 370, 20140124
  * Menéndez-Velázquez et al. (2022) Menéndez-Velázquez A., Morales D., García-Delgado A. B., 2022, Int J Environ Res Public Health, 19
  * Noll et al. (2012) Noll S., Kausch W., Barden M., Jones A. M., Szyszka C., Kimeswenger S., Vinther J., 2012, A&A, 543, A92
  * Owens et al. (2020) Owens A. C., Cochard P., Durrant J., Farnworth B., Perkin E. K., Seymoure B., 2020, Biological Conservation, 241, 108259
  * Parkinson et al. (2020) Parkinson E., Lawson J., Tiegs S. D., 2020, PLOS ONE, 15, 1
  * Patat et al. (2006) Patat F., Ugolnikov O. S., Postylyakov O. V., 2006, Astronomy and Astrophysics, 455, 385
  * Perkin et al. (2011) Perkin E. K., Hölker F., Richardson J. S., Sadler J. P., Wolter C., Tockner K., 2011, Ecosphere, 2, art122
  * Perkin et al. (2014) Perkin E. K., Hölker F., Tockner K., 2014, Freshwater Biology, 59, 368
  * Posch et al. (2018) Posch T., Binder F., Puschnig J., 2018, Journal of Quantitative Spectroscopy and Radiative Transfer, 211, 144
  * Puschnig et al. (2014a) Puschnig J., Posch T., Uttenthaler S., 2014a, Journal of Quantitative Spectroscopy and Radiative Transfer, 139, 64
  * Puschnig et al. (2014b) Puschnig J., Schwope A., Posch T., Schwarz R., 2014b, Journal of Quantitative Spectroscopy and Radiative Transfer, 139, 76
  * Puschnig et al. (2020) Puschnig J., Wallner S., Posch T., 2020, MNRAS, 492, 2622
  * Puschnig et al. (2021) Puschnig J., Näslund M., Schwope A., Wallner S., 2021, Monthly Notices of the Royal Astronomical Society, 502, 1095
  * Roesch et al. (2002) Roesch A., Wild M., Pinker R., Ohmura A., 2002, Journal of Geophysical Research: Atmospheres, 107, ACL 13
  * Schmidt & Spoelstra (2020) Schmidt W., Spoelstra H., 2020, Darkness monitoring in the Netherlands 2009-2019, doi:10.5281/zenodo.4293366, https://doi.org/10.5281/zenodo.4293366
  * Schwarz (1978) Schwarz G., 1978, The Annals of Statistics, 6, 461
  * Schwarz et al. (2006) Schwarz J. P., et al., 2006, Journal of Geophysical Research: Atmospheres, 111
  * Ściężor (2020) Ściężor T., 2020, Journal of Quantitative Spectroscopy and Radiative Transfer, 247, 106962
  * Ściężor & Czaplicka (2020) Ściężor T., Czaplicka A., 2020, Journal of Quantitative Spectroscopy and Radiative Transfer, 254, 107168
  * Ściężor & Kubala (2014) Ściężor T., Kubala M., 2014, Monthly Notices of the Royal Astronomical Society, 444, 2487
  * Tibshirani (1996) Tibshirani R., 1996, Journal of the Royal Statistical Society: Series B (Methodological), 58, 267
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science & Engineering, 13, 22
  * Wallner & Kocifaj (2019) Wallner S., Kocifaj M., 2019, J. Quant. Spectrosc. Radiative Transfer, 239, 106648
  * Zamorano et al. (2015) Zamorano J., et al., 2015, IAU General Assembly, 22, 2254642

##  Appendix A A synthetic sky model for the SQM band as a Reference

In order to quantify anthropogenic light at night, knowledge of the natural
night sky brightness is needed. An all-sky model that takes into account
scattered moonlight, starlight, molecular emission of the lower atmosphere,
emission lines in the upper atmosphere and airglow, was published by Noll et
al. (2012) and Jones et al. (2013), as part of an Austrian in-kind
contribution to the European Southern Observatory (ESO), e.g. ESO’s exposure
time calculator is based on the model.

The main input parameters are zenith distance or airmass of the observation,
precipitable water vapor and monthly averaged solar flux. For the moon
radiance component, the separation of Sun and Moon as seen from Earth, the
Moon-target separation, Moon altitude over horizon and the Moon-Earth distance
are needed.

The result is a synthetic (cloud-free) night sky spectrum for the target
location.

For our purpose, we multiply the so derived night sky spectrum with the SQM
transmission curve and calculate a synthetic SQM magnitude via integration. We
have decided to make some simplifications, allowing us to evaluate the model
on a 2-dimensional parameter grid with vectors of (Sun-Moon-separation, Moon
altitude) only. This is reasonable in our case, because the measurement
devices we are using, the Sky Quality Meters of type SQM-LE, are equipped with
a front lens that narrows down the field of view to a roughly 20 degree wide
cone, pointed towards zenith. Hence, we only need to consider zenithal night
sky brightness. The two input parameters Moon-target separation and Moon
altitude can thus be simplified to one parameter, with the former one being
the Moon zenith distance. We have further decided to evaluate the model for a
fixed precipitable water vapor value of 5 mm, a monthly averaged solar flux of
130 sfu and for a fixed mean Moon-Earth distance. These simplifications have
practically no influence on our results, since ALAN’s contribution to our SQM
measurements is magnitudes larger than the natural variation caused by
phenomena such as water vapor or solar flux. However, variations due to Moon
phase and height are large enough to be important and are thus fully treated
by our gridded model evaluation for the zenith.

Since the natural, cloudless sky brightness changes smoothly, a grid spacing
of one degree in both parameters (Sun-Moon-separation, Moon altitude) was
found to be sufficient.

The Sky Quality Meters (SQMs) are produced by
Unihedron444http://www.unihedron.com/. Several types are available, of which
we are using those denoted by SQM-LE, indicating that they are equipped with a
lens (L), and connected via ethernet (E). The lens narrows down the field of
view to a half-width-half-maximum (HWHM) angular sensitivity of 10 degree,
corresponding to a cone with an opening angle of 20 degree. The spectral
sensitivity of the SQM is sometimes compared to the V-band, but as shown in
Figure 6 it is much wider towards the blue end. Detailed technical
specifications are found in Cinzano (2005).

All our SQMs are operated inside weather-proof housings, also produced by
Unihedron. Our measurements are corrected for the loss of light due to the
housing’s cover glass (via subtraction of 0.11
${\mathrm{mag}_{\mathrm{SQM}}\,\mathrm{arcsec}^{-2}}$ from the reading).

We aim to compare our SQM measurements to natural (cloud-free) night sky
brightness values derived through a sky model. From our gridded model
implementation, spectra in physical flux units of $\mathrm{photons}\ s^{-1}\
m^{-2}\ \mathrm{\upmu m}^{-1}\ \mathrm{arcsec}^{-2}$ are obtained, which we
convert to $\mathrm{erg}\ s^{-1}\ \mathrm{cm}^{-2}\ \text{\AA}^{-1}\
\mathrm{arcsec}^{-2}$. The spectra are then multiplied with the SQM
transmission curve. Subsequently, the total flux is calculated via integration
and the result is divided by the bandwidth, which finally gives again a flux
density in $\mathrm{erg}\ s^{-1}\ \mathrm{cm}^{-2}\ \text{\AA}^{-1}\
\mathrm{arcsec}^{-2}$. The transformation into a Vega-based SQM magnitude
(that should be compared to our measurements), is done using the photometric
zeropoint and the zero magnitude flux for the SQM. The latter two quantities
were derived as explained in the following:

  * •

A Vega spectrum was downloaded from the CALSPEC Calibration
Database555https://archive.stsci.edu/hlsps/reference-
atlases/cdbs/current_calspec/. The flux density of the spectrum is given in
units of $\mathrm{erg}\ s^{-1}\ \mathrm{cm}^{-2}\ \text{\AA}^{-1}$. For Vega
(A0V star) the V-band magnitude is 0.03 mag and its (B-V) color is 0.

  * •

The Vega spectrum is multiplied with the SQM transmission curve, and for
sanity checks also with B, V, R Bessel filter transmission curves using data
from Bessell (1990).

  * •

The bandwidths were calculated via integration of the filter curves along the
wavelength axis (see Figure 6 for the transmission curves used). As seen in
Table 4, our results agree on a 10-percent level with those published by
Masana et al. (2022), who found 959, 909, 1634 and 2228 Å for B, V, R and the
SQM band.

  * •

The integrated flux of the product (Vega*Filter) is divided by the bandwidth.
The resulting flux density (FD) is then given in units of $\mathrm{erg}\
s^{-1}\ \mathrm{cm}^{-2}\ \text{\AA}^{-1}$. For B, V, R the zeropoints (ZP)
are calculated via: $0.03=-2.5\mathrm{log(FD)+ZP}$. Bessell (2005) published
zeropoints for these filters. They are -20.45, -21.12 and -21.61 mag for B, V
and R respectively. Our results agree at a level of 0.1 mag, which is better
than the absolute photometric accuracy of the SQM.

  * •

It is known that Unihedron’s in-house calibration of the SQM is performed
using a green LED. This calibration leads to (SQM$-$V) = $-$0.35 (see Cinzano,
2005, Figure 15). Thus, the SQM photometric zeropoint is calculated via:
$-0.32=-2.5\mathrm{log(FD)+ZP}$. The SQM zeropoint is -21.18$\pm$0.1 mag. The
zeropoint is defined as the magnitude of a source having unity flux.

  * •

Finally, the zero magnitude flux (FD0) can be calculated using the zeropoint
from above and: $0=-2.5\mathrm{log(FD_{0})+ZP}$. For the SQM the zero
magnitude flux is $3.38\ \times\ 10^{-9}$ $\mathrm{erg}\ s^{-1}\
\mathrm{cm}^{-2}\ \text{\AA}^{-1}$.

Table 4: Filter (column 1), peak wavelength (column 2) and bandwidth (column 3). The zeropoint (ZP), defined as the magnitude of a source having a flux of 1 $\mathrm{erg}\ s^{-1}\ \mathrm{cm}^{-2}\ \text{\AA}^{-1}$(column 4), and the zero-magnitude flux (column 5), were derived via multiplication with subsequent integration of a Vega spectrum with the corresponding filter transmission curve. Filter | $\lambda_{peak}$ | BW | ZP | $FD_{0}$
---|---|---|---|---
| [Å] | [Å] | [$\mathrm{mag}$] | [$10^{-9}$ $\mathrm{erg}\ s^{-1}\ \mathrm{cm}^{-2}\ \text{\AA}^{-1}$]
B | 4350 | 959 | -20.48 | 6.43
V | 5437 | 893 | -21.08 | 3.69
R | 6433 | 1591 | -21.64 | 2.21
SQM | 5087 | 2008 | -21.18 | 3.38
Figure 6: Bessel B, V, R and SQM transmission curves (solid curves) and
Gaussian fit results (dashed curves) on top of the spectrum of Vega.

##  Appendix B Long-term trends

Figure 7: Long-term trend for STO and IFA. See caption Figure 3 for more
details.

Figure 8: Long-term trend for GRA and FOA. See caption Figure 3 for more
details.

Figure 9: Long-term trend for BRA and FEU. See caption Figure 3 for more
details.

Figure 10: Long-term trend for FRE and GIS. See caption Figure 3 for more
details.

Figure 11: Long-term trend for GRI and GRU. See caption Figure 3 for more
details.

Figure 12: Long-term trend for KID and KRI. See caption Figure 3 for more
details.

Figure 13: Long-term trend for LGO and LSM. See caption Figure 3 for more
details.

Figure 14: Long-term trend for LOS and MAT. See caption Figure 3 for more
details.

Figure 15: Long-term trend for MUN and PAS. See caption Figure 3 for more
details.

Figure 16: Long-term trend for STY and TRA. See caption Figure 3 for more
details.

Figure 17: Long-term trend for ULI and VOE. See caption Figure 3 for more
details.

Figure 18: Long-term trend for WEL and ZOE. See caption Figure 3 for more
details.