Analysis of Outdoor and Indoor Radon Concentration Time Series Recorded with RadonEye Monitors

1. Introduction

1.1. Radon as Pollutant and as Tracer

Radon (Rn) is widely known as hazardous pollutant, which can cause lung cancer [1]. Consequently, there have been great efforts globally to monitor indoor Rn concentrations, assess levels of risk, and optimize Rn abatement policy and practice. As a direct result of these efforts there is now extensive national and international regulation of indoor Rn, e.g., the European Basic Safety Standards [2]. In recent decades the ever-present public health risk posed by Rn has resulted in the development of a wide variety of increasingly affordable and portable passive and active Rn monitors, e.g. [3,4,5,6,7,8,9,10,11,12].

Consequently, large databases of indoor and outdoor Rn measurements have been developed using these instruments, but (i) this data is not always easily available due to data protection issues, and (ii) the sensitivity and/or temporal resolution of these measurements may not be particularly high because reference levels for Rn are typically between 100 – 300 Bq m${}^{-3}$ as annual mean.

But Rn also has a second face. Due to its unique collection of physical properties, the noble gas Rn can also serve as a convenient, efficient tracer in environmental studies of transport and mixing processes in the ground, hydrosphere and atmosphere [13,14,15]. One main technique in tracer research is the analysis of Rn concentration time series. Such studies seek to identify "signals" associated with some natural events and attempt to relate these signals to observed Rn variability. For atmospheric studies in particular, environmental Rn concentrations can be orders of magnitude below indoor or sub-surface radon concentrations, and exhibit variability on much shorter timescales (minutes to hours). These needs have driven the development of a separate collection of "direct" and "by progeny" active Rn detectors capable of measuring concentrations below 1 Bq m${}^{-3}$ at sub-hourly temporal resolution (e.g., [16,17,18,19,20,21,22]). Large databases of outdoor Rn observations are now also being accumulated using these instruments associated with the WMO GAW and ICOS research networks. Unfortunately, it is not always possible to utilise outdoor Rn datasets collected by public health Rn monitors in the same ways as research grade environmental Rn monitors due to differences in sensitivity, precision, and temporal resolution. However, due to the much greater affordability of public health Rn monitors, there is value in increasing our understanding of their measurement performance in all settings.

1.2. Radon Source, Transport, and Temporal Variability

Environmental radon (Rn; here limited to the isotope ${}^{222}$Rn from the ${}^{238}$U decay series) has a predominantly terrestrial source. Rn emissions from ocean surfaces are typically 2-to-3 orders of magnitude less [23,24]. Its parent nuclide, ${}^{226}$Ra, is ubiquitous in soils, rocks, and most building materials. While the formation of Rn from Ra is relatively constant in time ($T_{1/2}$${}^{226}$Ra ∼1600 y), the release of Rn from soil grains, as well as its transport within and from the ground, is strongly influenced by changing soil moisture conditions and soil freezing [25,26,27]. Once Rn has entered the atmosphere, it is readily advected by wind. Because of the global Rn source heterogeneity, Rn’s 3.8 d half-life, and diurnally changing atmospheric mixing conditions, environmental Rn concentrations can show strong temporal variability which, in general, is difficult to model precisely. The reasons for this difficulty include missing - unavailable or unknown – predictors, and the overall high complexity of the "environmental Rn" system (see, e.g., the scheme "From rock to risk", [13]).

While Rn enters the atmosphere predominantly by diffusive transport, once in the atmosphere it is distributed by a combination of advective, turbulent and diffusive transport, depending on prevailing atmospheric conditions. In the absence of active convective clouds or strong fronts ([28,29]), Rn is essentially confined to the planetary boundary layer (PBL), the vertical extent of which is controlled at any given time by the combined total of convective and mechanical mixing processes. At night or early morning over land, when convective mixing is suppressed, the PBL depth can be quite shallow if wind speeds are light, resulting in higher Rn concentrations near the ground. When wind speeds increase, or daytime convection is active, the Rn concentration near the ground decreases due to dilution. The combined contributions of mechanical and convective turbulence at any given time dictate the atmospheric mixing state. One popular measure of the atmospheric mixing state (or stability) is the Richardson number, which is the ratio between buoyancy influences (driving convection) and the vertical wind speed gradient (driving mechanical mixing, in combination with the surface roughness). Diffusion of Rn in the atmosphere plays little role unless the atmosphere is very stable.

1.3. Low-Cost Radon Measurement and Citizen Science

Consumer grade radon monitors, which are orders of magnitude cheaper than professional "research grade" instruments, have become increasingly popular. To date, they have mainly been used in the context of private screening measurements and for Citizen Science. At the time of writing this article, the most sensitive among the consumer-grade radon monitors was the RadonEye monitor (FT Lab, South Korea [30]). Due to its economical price (less than 200€ for the basic version) it has become very popular and has been widely adopted by private and professional users, including use in Citizen Science projects. The RadonEye has been subjected to a number of QA studies (see section 2.1), the results of which appear favorable. Nevertheless, it is important to understand its limitations and define the objectives for its use. To this end, it is the intention of this study to scrutinize the performance of this monitor for a selection of indoor and outdoor applications.

1.4. Objective of the Paper

The objectives of this work are (1) to improve understanding of certain aspects of the counting statistics for results acquired by RadonEye monitors, and (2) to perform descriptive and exploratory statistics of indoor and outdoor Rn concentration timeseries collected using a RadonEye.

(1) The RadonEye only reports Rn concentrations in Bq m${}^{-3}$ rounded to integers, leading to questions about the expected true concentration, conditional to the reported value. Since the exact algorithm used by the RadonEye is not publicly-available, we seek to better understand the reported results. Furthermore, we estimate detection limits and confidence intervals of the RadonEye count rates.

(2) We analyze indoor and outdoor Rn concentration timeseries recorded by RadonEye monitors over periods of several months. In addition to basic descriptive statistics, we apply autoregressive and fractal analysis, which has so far been less common in environmental radon studies. Tools are traditional time-series decomposition (analysis of trend and periodicity) and methods from the AR family (auto-regressive analysis). Further insight may be gained using multifractal methods such as the Hurst coefficient and $f\left(\alpha \right)$ spectra. The latter may help recovering long-term memory structures, which from a physics perspective reflect meteorological episodes (weather periods) rather than short-term variability.

The paper is organized as follows: the methods section has been split into two parts. The first part deals with measurement hardware and locations. The second part is about statistical procedures. This is followed by the results, with two sub-sections related to (1) and (2), above. Discussion of the findings and suggestions for further research are presented in the final section.

2. Methods 1: Hardware, Experimental Design and Data

2.1. The RadonEye Monitor

The instrument used to generate Rn time series is the RadonEye, based on an ionization chamber with sensitivity 1.35 cpm/(100 Bq m${}^{-3}$) = 0.81 cph/(Bq m${}^{-3}$) according to the manufacturer, or inversely, 1 cph corresponds to about 1.2 Bq m${}^{-3}$. It is not clear whether this is a generic value applied to all exemplars or a mean over individual calibration factors, and whether in each device the generic or an individual calibration is implemented. Its output are Rn concentrations which can be retrieved by smartphone app via Bluetooth. 10 min means are displayed but only 1 h-means rounded to integers are stored in a simple text file. QA studies - not to be discussed here in detail - show acceptable performance for usual indoor Rn concentrations [8,9] ; investigations with very low Rn concentrations (outdoors: typically <1 - 30 Bq m${}^{-3}$) are still missing. Problems seem to exist with sensitivity against ${}^{220}$Rn (thoron, Tn) [31] and dependence on humidity which may introduce bias, as well as with performance at extreme temperatures. Apparently for the latter reason the manufacturer does not recommend its use outdoors. Since the objective of this study is not the absolute value of Rn concentration but its dynamic, we shall not focus further on this possible problem. (However also the temporal behaviour may be impaired to some degree; this has however not yet been investigated.) One drawback is the lack of an internal buffer battery: in case of power loss the unsaved part of the recorded series is lost and a new series is started. Pictures of RadonEyes are shown in Figure 1.

Investigations concerning counting statistics and how to interpret reported results are presented in section 4.

Ionization chambers suffer from internal background (BG) due to the accumulation of the long-lived alpha emitting members of the ${}^{222}$Rn decay chain ${}^{210}$Po and ${}^{210}$Pb. Newly acquired devices, as used in this project, should therefore have very low internal BG. It seems that for new devices, one has to expect the equivalent of 1 Bq m${}^{-3}$ or less (personal communication, BfS Berlin). This may be compared with sensitivity, 0.81 cph/(Bq m${}^{-3}$) or less than 1 count per hour (which is the reporting period) for exposure to 1 Bq m${}^{-3}$, which is common for outdoor Rn.

2.2. Measurement Locations and Periods

Vienna: Indoor radon concentration was measured hourly in Vienna, Austria, between 17 February and 23 July 2024. The monitor was placed on a rack near a bookshelf about 50 cm above the wooden floor. The room is located on the second floor of an apartment building built in the late 19th or early 20th century, made of brick, typical for Vienna. The building has four floors plus a basement whose floor is not isolated against the ground and which is also not isolated from the interior of the building. The room is about 3.6 m high and has an area of about 20 m${}^2$. Two windows look onto a courtyard; their wooden frames are quite old and, in spite of double glass, they are not very airtight. The local geology is Holocene alluvial of the Danube [32] with low geogenic radon potential. Previous measurements in the same room (19 August to 12 November 2000) using a Pylon-AB4 monitor (based on Lucas cell) consisting of 2093 hourly means (as for the present analysis) yielded a mean Rn concentration of 15 Bq m${}^{-3}$ only (min-max = 0 - 46). However, windows were open during part of the measurement period. If the windows were closed, the concentration rose to about twice as high, indicating influence of either geogenic Rn or exhalation from building material. Another 2-week trial measurement by RadonEye in April 2023 gave 30 Bq m${}^{-3}$ (windows closed). The location is 48.22° N latitude / 16.37° E longitude.

Berlin: Outdoor Rn was measured in Berlin, Germany between 23 February and 29 July 2024. The monitor was placed in a protective plastic case Figure (Figure 1) on a wooden balcony about 2 m away from the wall and about 6 m above ground, a large garden. The ground is post-Weichsel glacial alluvial of the Spree Urstromtal [33], mainly sand, with low geogenic radon potential. Long-term measurement using track-etch passive detectors at a nearby location, 1.5 m above ground, gave a 3-year-mean outdoor concentration of 6 Bq m${}^{-3}$[34]. Four test measurement series with RadonEye in May, July/August, September and October 2023 gave 3, 3, 7 and 6 Bq m${}^{-3}$, respectively.

Berlin: Indoor Rn: Two indoor Rn series were recorded in Berlin. The first between 25 October 2023 and 3 February 2024. In the same apartment, radon has been measured indoors partly in parallel with outdoor Rn between 2 June and 29 July 2024. The RadonEye was placed like shown in Figure 1, left, during the first (winter) period and in a plastic box open upwards, located nearby directly on the tile floor, during the second (summer) period. Possible thoron influence is probably negligible. During measurement, windows were mostly closed. The apartment is located below the wooden roof (corresponding second floor) of a building about 90-100 years old. There is a souterrain-type basement with concrete floor, normally separated from the interior of the house by a locked, but not air tight door. The location is 52.49 N latitude / 13.53 E longitude.

3. Methods 2: Statistical Analysis of Reported Rn Time Series

Basic methods (descriptive and exploratory statistics, periodogram) follow established rules and will not be discussed here. This section is concerned with the interpretation of the reported results, not with measurement statistics. This is discussed in Section 4.

3.1. Box-Jenkins Scheme

3.1.1. Autoregressive Moving Average Model

This procedure attempts to separate variability components of a time series, test for the presence of a trend, adjust the series by removing trend and periodicity thus identify aperiodic variability. Further it helps estimating an autoregressive model that allows prediction. Thus, the residuals about the model can be traced to the aperiodic components. The most comprehensive textbook is [35]. For a short introduction, e.g., [36], the NIST handbook [37], section 6.4.4.5 and the Wikipedia entry "Box–Jenkins method". Time series are conceived as realizations of autoregressive (AR) processes,

$$x_t=m+{\alpha }_1{x}_{t-1}+{\alpha }_2{x}_{t-2}+...$$

(1)

or moving average (MA) processes,

$$x_t=m+{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2}+...$$

(2)

where $\epsilon$ are $i.i.d.$ random terms associated to the i-th member of the series; thus the value at t is influenced by the noise of pdimension revious terms. The highest index of $\alpha$ and $\theta$ is the order of the process. AR and AM are combined to ARMA($p,q$) processes, p and q the AR and MA orders. The concept is extended to ARIMA($p,d,q$) or autoregressive integrated moving average process. The order d can be fractional (F), resulting in ARFIMA or FARIMA [38] models. The parameter d refers to differencing the original series, $y_t=x_t-x_{t-1}$, repeated d times (fractional d require particular definition), and applying the ARMA procedure to Y. For application of the ARFIMA technique to Rn series, see [39]. Differencing X is performed if the original series is non-stationary. The technique is much used in financial mathematics, in climatology and other disciplines which deal with time series, but less so in the field of Rn. As one notable paper we mention Stransky [40], who applied ARIMA with external (exogenous) predictors to a Rn time series recorded in a cave.

3.1.2. Box-Jenkins Scheme - Autocorrelation and Partial Autocorrelation Functions (ACF and PACF)

Analysis of the ACF and PACF enables estimation of the model structure in terms of AR and AM order. The PACF removes the influence of correlations for each preceding lag., see [35], section 3.2.5. How this is done to estimate the orders of ARMA($p,q$) will be shown in the results Section 5.2.2.

The correlation length, that is the period during which the values are considered correlated, can be quantified by the time after which the correlation decreases below a threshold. Commonly used ones are $1/e=0.37$ or $p=0.05$ as error probability of the null hypothesis: "not correlated".

3.2. Stationarity

Stationarity is a rather tricky concept and there are different definitions and approaches. A stochastic process X is called stationary if $E\left[X\right]=const$. Another definition is $E[X+t]=E\left[X\right]$ for all t, that is, time shift does not change the expectation. In general this cannot be proved because usually one has only one realization $x_1,...,x_n$ of X, so that the about equivalent Birkhof-ergodicity, loosely put $E_t\left[X_{\omega }\right]\forall \omega =E_{\omega }\left[X_t\right]\forall t$ cannot be tested ($\omega$ - index that counts realizations of process X). Finding a trend within a limited domain (the observation time or space) does not indicate non-stationarity because in a larger domain the picture may be different. For example, the indoor Rn series in Berlin both look about stationary, but they are different between winter and summer, indicating an annual cycle (which exists as is known from other experiments). But extending the series over many years would again let appear it stationary. Non-stationarity would be given if for some physical reason, such as climatic change, there would be a long-term increasing or falling trend; but such effect has so far not been identified. For the purpose of our analysis we may assume the series stationary. A non-decaying ACF indicates that the series is probably non-stationary in the investigated domain.

To complicate further, stationarity and stationarity on top of a trend has to be distinguished. In the latter case the presence of a deterministic trend is assumed, around which the series fluctuates. The question is whether this fluctuation is random walk-type, i.e. non stationary.

In spite of the assumption of stationarity based on underlying physics, some statistical tests were applied.

The non-parametric Mann-Kendall test [41] (see also the Past4 manual) looks for a deterministic trend in the series. The Dickey-Fuller stationarity test [42] is based on differenced series; it tests whether coefficients of differences are equal 1, therefore it belongs to so-called unit root tests. the p-value of the 0-hypothesis root=1 is determined, which indicates non-stationarity. Here the basic and the augmented versions of the test as implemented in Gretl [43] were used.

3.3. Multifractal Methods

Fractal and multifractal aim to explore the "irregular" behaviour of a field or a time series, that is the variability beyond trend and periodicity. This manifests in non-periodic "episodes" of higher or lower values of the response variable (here, Rn concentration) and in the "roughness" of the series graph.

3.3.1. Hurst Coefficient

The term Long Memory "refers to persisting correlation between distant observations in a time series" [44], excluding correlation due to periodicity. Long memory indicates the presence of aperiodic episodes in the time series, typical for Rn series, owing to weather episodes as meteorology is the most important predictor. If one succeeds to explain the episodes in the Rn series from meteorological episodes, the remaining unexplained variability indicates the presence of additional predictors not accounted for. Many authors attempted to use signals in ground and hydrological Rn time series for seismic prediction because it is empirically known that seismic / tectonic phenomena can produce Rn signals. However, no reliable seismic prediction based on Rn has yet been successful, if one understands acceptably low first and second kind prediction error probabilities as criterion of reliability.

The Hurst coefficient is a measure of persistence, or the tendency that the variable approximately maintains a certain level for some time. 0.5<H<1, with higher H meaning less volatile or smoother series. High (low) values tend to be followed by high (low) values, which can be interpreted by (irregularly) recurring "episodes" which last some time. H<0.5 indicates an anti-persistent process, meaning that values tend to change quickly between low and high. H=0.5 indicates Brownian motion. Bowers & Tung [38] cite a number of physical phenomena for which long-term dependence has been found. An interesting paper of the history of the Hurst effect is Graves [45], which also discusses the connection with ARFIMA.

For an overview of methods to estimate the Hurst coefficient, [46,47,48]. A critical review is given in Donner [49]; it is not clear to which degree it applies to the methods used in this paper.

3.3.2. Multifractal $f\left(\alpha \right)$ Spectrum

Multifractals can be understood as interwoven fractals of point sets which are realizations of random fields, where an "intensity" is given at locations of the carrier domain (here: the time axis). This can be understood as the density of a measure (in the mathematical sense) within a "box" (time interval) of given size h. It is assumed that locally, the field behaves as $\left|Y(t+h)-Y\left(t\right)\right|\sim h^{\alpha \left(t\right)}$, where $\alpha \left(t\right)$ is called the local Hölder exponent or singularity strength at location t. The number of boxes $n\left(\alpha \right)d\alpha$ with Hölder exponent $\in (\alpha ,\alpha +d\alpha )$behaves as $n\left(\alpha \right)\sim h^{-f\left(\alpha \right)}$, and the graph $(\alpha ,f(\alpha \left)\right)$ is called multifractal spectrum. It looks like an upside down parabola which is the wider, the "rougher" the graph of the time series is (The same applies to spatial fields). The dimension $f\left(\alpha \right)$ for high $\alpha$ characterizes the pattern of extremes of Y, in this case the Rn concentration. For rigid definitions and a more detailed discussion see, e.g., [50,51] or the textbook [52].

The spectra were calculated with the algorithm proposed by [50], which essentially rests on linear regression of log(partition functions) against log(scale). A first version of the multifractal analysis of spatial and temporal radon fields was presented in [53].

3.3.3. Fractal Dimension of the Graph

The graph interpreted as a point set has a fractal dimension D which measures its "roughness". Among infinitely many possible definitions of D usually the box dimension $D_B$ is chosen in this context. Essentially, the graph space $(x,f(x\left)\right)$, with embedding Euclidean dimension 2 in our case is divided into equal boxes and the number of boxes "occupied" by points of the graph is counted in dependence of box size or scale. The number scales like $N\sim {\delta }^{-D}$ where the scale $\delta$ is the box size relative to some initial size or as a fraction of the area of the entire domain. D can therefore be estimated as slope of the $logN/log\delta$ plot. By theory, $D=2-H$, with 2 = the embedding dimension of a time series graph.

For analyses related to Rn, see [47] who found good agreement between estimations of D via box counting and the Hurst coefficient.

3.3.4. Garger Exponent

Garger et al. [54] investigated the fluctuation of radionuclide concentrations in aerosols after the Chernobyl accident (1986). They found a power law dependence of the fluctuation over a period from this time, $s\left(T\right)\sim T^{-g}$, for not too long aggregation periods T. They found $g\approx 0.33$ which would coincide with what one expects for turbulent wind fields. The approach was further elaborated by Hatano et al. [55,56]. Smooth functions $X\left(T\right)$ would yield $g=1$, white noise, $g=0$ and random walk, $g=0.5$. Concerning Rn, the coefficient was first calculated in 2002 [57] for the several indoor Rn series, resulting quite consistently in g = 0.60 - 0.62.

$s\left(T\right)$ is defined, $s\left(T\right):=SD\left(X\right(T\left)\right)$ with $X\left(T\right)=mean$ over series of length T.

3.3.5. Attractor Embedding Dimension

The phase space of the process can be reconstructed from the time series with an algorithm based on Taken’s embedding theorem; see [58] sec.5 and for an easy explanation, e.g., [59,60]. From the time series $x\left(t_{i}i\right)$, i = 1...n, we construct the E-dimensional series $y\left(t_i\right):=(x\left(t_i\right),x(t_i+\tau ),...,x(t_i+(E-1)\tau ))$, where E is the Euclidean embedding dimension and $\tau$ a time delay which has to be chosen appropriately that $x\left(t_i\right)$ and $x(t_i+\tau )$ are little correlated. A reasonable choice of $\tau$ is the decorrelation time according to the ACF. The embedding dimension E can be seen as estimate of the number of variables which control the system.

Estimation of the minimal E is performed with two methods: (1) in the y-space, which is topologically equivalent to the phase space, the correlation dimension $D_c$ is determined: For every $y\left(t_i\right)$, one counts the number N of points within distance $ϵ$. The mean $A{M}_i\left(N\right)$ scales as $AM\left(N\right)\sim {ϵ}^{Dc}$, so that $D_c$ can be estimated as slope of the $log\left(AM\right)/logϵ$ plot. The correct E is found when $D_c$ levels off with increasing E. (2) In the false nearest neighbor (FNN) technique one chooses E such as to minimize the number of so-called false near neighbors, which are points of the y-series which falsely appear as neighbors of $y\left(i\right)$ due to E set too small. One increases E until the FNN rate is reasonably low.

The lowest E and the maximal Lyapunov exponent were estimated with CDA and with R, packages "tseriesChaos" [61] and "nonlinearTseries" [62]. For estimation methods see [63].

3.3.6. Lyapunov Exponent

A qualitative definition is given in [64]: "The Lyapunov exponent (...) is (...) the measure of the sensitivity of a system to its initial conditions. Lyapunov exponents (...) are related to the exponentially fast divergence or convergence of nearby orbits in phase space. A positive Lyapunov exponent is an indication of chaotic behaviour in the time series, while a negative value of Lyapunov exponent shows periodic trajectories, or non-chaotic behaviour, of a dynamical system. The more positive the value of Lyapunov exponent is, the shorter is the time interval up to which the system’s future can be predicted and vice versa." Nearby orbits or trajectories are followed with time and their separation is measured. The maximal Lyapunov exponent is the rate of an anticipated exponentially growing separation. In [57] estimation of Lyapunov exponents together with phase space reconstruction through the recurrence plot technique, introduced by [56] and much used since, have been attempted, but little further literature in the field of Rn studies exists, quoted in the results section.

Software: basic operations and graphs were produced with Excel; for time series analysis, Past v. 4.03 and 4.17 [65], Statgraphics v.4 and Gretl 2019a [43], CDA 1.1 [66] and different R packages were used, detailed in the respective results sections. The Hurst exponent, the $f\left(\alpha \right)$ spectra were additionally computed with in-house developed software written in QB64 [67].

4. Results 1: Aspects of Counting Statistics

Basic properties of the RadonEye monitor were presented in Section 2.1. In this section, we shall address some issues related to low count rates and to the (unknown) internal evaluation algorithm of the device.

4.1. Critical Limits

The simple Currie approach [68,69] gives $decisionlimit=Lc\equiv x^{*}=z_{\alpha }\sqrt{2·BG}$, $detectionlimit=Ld\equiv x^{\#}=z_{\alpha }^2+2k_{\alpha }\sqrt{2·BG}$. With usual $k_{0.05}$ 1-sided = 1.645, this is $x^{*}=2.3\sqrt{BG}$ and $x^{\#}=2.71+4.65\sqrt{BG}$ For the newer ISO 11929 version, which yields slightly different values see, e.g., [70,71,72]. Here we stick with the Currie version for simplicity.

Assuming a BG of 1 count per hour, BG=1, we find $x^{*}$ = 2.3 counts = 2.9 Bq m${}^{-3}$ and $x^{\#}$ = 7.4 counts = 9.1 Bq m${}^{-3}$. On the other hand, in many instances 0 or a few Bq m${}^{-3}$ are reported; even if BG=1 may be too high for new RadonEyes, this shows that very low reported values should be understood cum grano salis. Dependence of the critical limits on BG is shown in Figure 2.

4.2. Frequentist Poisson Confidence Interval of Counts

The device counts pulses and its uncertainty can therefore be estimated from Poisson statistics. Figure 3a shows the dependence of relative uncertainty on concentration, for three measurement periods, under the Poisson hypothesis and the sensitivity reported by the manufacturer. Confidence intervals can be calculated by different methods, summarized in [73]. We prefer the simple ${\chi }_{\alpha /2;2x}^2/2,{\chi }_{(1-\alpha )/2;2(x+1)}^2/2$. Here ${\chi }^2$ is the inverse Chi-square with $2x$ and $2(x+1)$ degrees of freedom, respectively, x the count number. The $\alpha =0.05$ (i.e., 95%) - confidence limits for low count numbers are shown in Figure 3b. For x=0, CI95%=(0, 3.7), for x=1:(0.025, 5.6), x=2:(0.24, 7.2), and so on.

The often used Wald interval $(x+z_{\alpha /2}\sqrt{x},x+z_{1-\alpha /2}\sqrt{x})$, based on the normal approximation of the Poisson distribution, is applicable only for high observed counts x, say x>15. The reasoning would become more complicated if considering the BG.

In its smartphone App, the RadonEye reports a 95% confidence interval for the last one hour measurement based on six 10-min counting periods. The algorithm is not publicly available. For nominal 0 Bq m${}^{-3}$, a CI = (0, 19) is given, which appears too large to me if compared to the above estimate.

4.3. Bayesian Estimate of the Expected True Count Rate

We assume as prior the conjugate distribution of the Poisson distribution, which is the gamma distribution, $f\left(y\right)=[{\beta }^{\alpha }/\Gamma \left(\alpha \right)]y^{\alpha -1}{e}^{-\beta \gamma }$ (see, e.g., [74] or the Wikipedia entry "Poisson distribution". The predictive posterior distribution (i.e. the distribution of x) is the negative binomial. The posterior Poisson mean $m_{Bay}=(\alpha +x)/(\beta +n)$ for n=1, the number of observed counts of size x. We guess $\alpha$ and $\beta$ such that the gamma (i.e. the prior) mean $\alpha /\beta$ equals the midpoint of the frequentist confidence interval. Since Poisson variance = mean we find $\beta$=1 and $\alpha$=midpoint; however Bayesian inference has it that this choice is to an extent deliberate and may be understood as reflecting plausibility. $m_{Bay}$ is shown in Figure 3b. With $x\to \infty$ the Bayesian mean approaches the observed count number x. For x=0 observed counts, the estimated mean equals 0.92, for 1 count, m=1.9, and so on.

4.4. Bayesian Inversion of Reported Rounded Concentration Values

The RadonEye reports Rn concentrations rounded to integers, $C=cint(k/(calT\left)\right)$, $cal$ - calibration factor (cph/(Bq m${}^{-3}$)), k - counts per period T, $cint\left(y\right):=int(y+0.5)\equiv \lfloor y+0.5\rfloor \equiv floor(y+0.5)$ in different notation; therefore, re-converting into counts $k^{\prime }$ introduces additional uncertainty, $k^{\prime }=cint\left(Ccal\right)$. The actual rounding algorithm used in the RadonEye is not publicly available (For further information on the various rounding algorithms, please refer to e.g. the Wikipedia entry on "rounding").

The question is therefore, given a concentration C (rounded to integer) reported by the instrument, which is the expected true concentration or its most probable value (or more generally, its distribution)?

Let $C{\left(true\right)}_T$ the true mean concentration in interval T (here we focus on T=1 hour. In reality, Rn concentration is not entirely constant during one hour, see Figure 7, but considering this would render the reasoning very complicated. This may be left to a future study). In the following, the results refer to interval T. The index T is therefore suppressed.

C(true) leads to a number of pulses k within interval T. These are Poisson distributed with mean $\alpha =C·cal·T$, that is, $k\approx Poiss\left(\alpha \right)$. The probability to count k pulses is $p\left(k\right)=Poiss(\alpha ,k)={\alpha }^k{e}^{-\alpha }/k!.$

The registered concentration $C\left(reg\right)=k/(cal·T)$ and reported concentration, $C\left(rep\right)=cint\left(C\right(reg\left)\right)$. Some example values of $p\left(k\right)$ for true concentrations between 0.5 and 4.5 Bq m${}^{-3}$ are shown in Table 1. For instance, if the true concentration is 2.5 Bq m${}^{-3}$ (first row), the RadonEye reports 6 Bq m${}^{-3}$ (third column) with probability 3.7%. This table will be inverted in the next step.

$p\left(k\right)$ can be interpreted as conditional probability $p\left(C\right(rep\left)\right|C\left(true\right))$. The wanted posterior probability $p\left(C\right(true\left)\right|C\left(rep\right))$, which is the probability that the true concentration is $C\left(true\right)$ if $C\left(rep\right)$ - an integer number - is reported by the instrument, is therefore, by Bayes Law:

$p\left(C\left(true\right)\right|C\left(rep\right))=p\left(C\left(rep\right)\right|C\left(true\right))p\left(C\left(true\right)\right)/{\sum }_{C\left(true\right)}\left[p\left(C\left(rep\right)\right|C\left(true\right))p\left(C\left(true\right)\right)\right]$

For the prior $p\left(C\right(true\left)\right)$ we can choose a constant as an improper uninformative prior, so that it cancels in the numerator and denominator, and the desired posterior can readily be computed. Among many other possibilities, a sufficiently wide rectangular function centered around $C\left(rep\right)$ could also be chosen as a prior. (Of course, the estimated posterior always depends on the choice of prior to some extent.) From the posterior distribution, the conditional mean, the most probable value of the true concentration, and other statistics can be derived. Here, the calculation was done with a custom QBasic program.

The result using the improper-uninformative prior is shown in Figure 4. The most probable true concentration is practically the same as the reported. The expected true concentration is somewhat higher but the relative distance decreasing (lower right inset).

Re-converting the expected true concentration into counts and comparing with observed counts shows, like in the previous section, a decreasing difference (upper inset of Figure 4). Comparison with the more formal derivation in Section 4.3 shows very little difference. For observed x=0 and x=1 counts, the expected values equal 1.06 and 2.04, respectively, to be compared with 0.92 and 1.9.

4.5. Missing Nominal Concentration Values

Rounding induces a rather strange effect, namely, that certain integer values of $C\left(reg\right)$ can never occur, see Figure 5, for two values of cal and also in Table 1, third column. For cal=0.81 (default for the RadonEye), $C\left(rep\right)$ = 3, 8, 13,... Bq m${}^{-3}$ cannot occur. For cal=0.79, this would be $C\left(rep\right)$=2, 7, 12,.. Empirically, for three investigated RadonEyes missing nominal concentration values given in Table 2 are observed. We can see that 1) the missing values are different between the devices and 2) that the sequences are not regular as they should be if the proposed explanation is sufficient. Evidently it is not, although it is probably the core of the problem. A satisfactory explanation should be given by the manufacturer.

As a summary of the preceding sections, several open problems remain to be solved before the findings can be implemented in practice. However, only the manufacturer can provide detailed information about the algorithm. Therefore, the subsequent analyses of the empirical time series will rely on the reported raw data. The error is likely not significant, as the Bayes-expected concentrations, conditioned on the raw data, do not differ substantially from the raw values, except for concentrations close to zero. Moreover, this paper focuses more on the dynamic patterns of the series rather than absolute values, so a potential bias in these values should not significantly affect the results.

5. Results 2: Descriptive and Exploratory Statistics

5.1. Descriptive Statistics, Periodicity And Autocorrelation

5.1.1. Time Series and Descriptive Statistics

Table 3 presents basic statistics of the radon series graphically shown in Figure 6. The maxima observed in Berlin (indoor) about end January 2024 and Vienna end March / early April 2024 coincide with Sahara dust episodes. To our knowledge, the exact physical reason has not yet been explored. We hypothesize that Southern wind implies that air travels longer distances above territories with higher geogenic Rn exhalation than Westerly winds (predominant in Vienna and Berlin) that origin over the Atlantic Ocean. The short 10-day outdoor series recorded during a very hot period in mid-August 2024, Figure 7, shows the high diurnal variability, and that the maxima occur during morning hours. Low concentration values should be understood subject to Section 4.3 and Section 4.4.

5.1.2. Stationarity

The Mann-Kendall test for trend (performed with Past 4.04 because the more recent version 4.17 seems to have an issue with this test) indicates presence of a significant trend for all series. This may however be an accidental effect particular to the chosen measurement periods. According to the simple as well to the (stronger) augmented Dickey-Fuller test performed with Gretl, all series are approximately stationary at p<0.05 if the presence of a trend is considered.

5.1.3. Periodicity and Autocorrelation

Periodicity and autocorrelation (the ACF, autocorrelation function; the PACF - partial autocorrelation function - is discussed in Section 5.2) are summarized in Table 4 and visualized in Figure 8. The periodograms show the expected 24 h periodicity. In addition, one can recognize periodicities whose physical origin is unclear. The one of 26 - 27 days may point to lunar tides, although the synodic lunar period (time between the same phases of the moon) equals 29.53 days, the sidereal period (same position relative to the stars), 27.32 days. The other periodicities may reflect meteorological patterns: Periods of persistent inversion in winter and hot periods in summer with little wind and low turbulence are common in Central Europe (see, e.g., the several week long heat periods in summer 2024, interrupted only by rainy and stormy days which last few days, and the likewise week long notoriously depressing grey and cold inversion periods in winter). The diffuse clusters of periodicities around 3 days and 10 days may refer to such patterns. Similar interpretations were reported by Kikaj et al. [75]. The slope $\alpha$ of the log-log periodogram for small frequencies refers to the model: $power=f^{-\alpha }$, indicating the presence of coloured noise, see Section 5.3. For high frequencies, the periodograms show a flattening of the power with frequency, recognizable from the graphs in frequency domain (inset in the graph Figure 8, upper). In this range (temporal resolution of 10 h and below) the periodogram indicates that the reported data tends progressively to white noise. The probable reason is instrumental noise. This has also been found in [57] and in the TraceRadon project [76].

The autocorrelation plots show that in Vienna Rn reacts more slowly to external drivers than in Berlin indoors during the summer period, but similarly to Berlin indoors in winter (series Berlin-in 1). For some reason yet to be explained, indoor Rn dynamic in Berlin is different between the winter and summer periods: in winter the daily period is much less pronounced than in summer. Mean Rn concentration is higher in winter than in summer - a typical effect, perhaps due the stack effect induced by heating. This may also explain why it reacts to outdoor conditions less than in summer. Autocorrelation decreases very fast from lag 0 to lag 1, which "nugget effect" is due to the high statistical uncertainty of the values.

5.1.4. Cross-Correlation Indoor-Outdoor Radon

Cross-correlation analysis between the outdoor and indoor (summer) series shows that in Berlin, outdoor Rn is 6 hours ahead of indoor Rn (series 2) in average (Figure 9). The low overall correlation is due to the very high uncertainty of low Rn concentration values below decision and detection limit, see Section 2.1. Nevertheless, the Fourier and autocorrelation analyses indicate that the dynamic pattern of Rn concentration is adequately represented, despite the potential for uncertainty regarding the reported low values. It is important to note that the lag of 6 hours (lag=6h) represents an average value, which may not be constant over the entire period.

The ratio between outdoor and indoor Rn concentration in Berlin during the period when the monitors were working in parallel is shown in Figure 10. The 6-h lag has been considered. The plotted values are the ratios between means in 24-h sliding windows. The anomalous peak towards the end of July is not explained. The apartment was not occupied and windows were closed during the entire period. The mean ratio is equal to 0.39, which can be interpreted that outdoor Rn contributed about 39% to indoor Rn during that period. The remaining is owed to geogenic Rn and in a minor part to the (probably very little) exhalation from building material. Exhalation from water can be excluded during that period.

5.1.5. Outdoor Radon Diurnal Pattern

Accounting for the main 24-hour periodicity, outdoor concentration values were aggregated per hour of the day in Figure 11. One can clearly recognize the nocturnal radon accumulation effect which leads to maxima at 5 in the morning and minima around middle afternoon, in average. The physical reason is the higher atmospheric stability at night which leads to lower atmospheric mixing layer height and resulting higher Rn concentration. Similar graphs are shown in [77,78,79,80,81], or [76] Figure 15, p.44 for the EEC. An extensive study of these patterns is presesnted in [82]. The time of the maximum and the onset of increasing and decreasing depends on sunshine hours. The series was recorded between February and July in Berlin, where sunrise varies between 7:05 and 4:20 CET in the morning during the measurement period, so that the curve is "smeared" as reflecting this variability. This is corroborated by Figure 12, which shows the same curve for the measurement months: one can notice that in tendency, the maxima and onsets shift to earlier and later hours, respectively, with increasing month. Kikaj et al. [83] Figure2, have shown that the actual shape of the diurnal Rn concentration curve depends on atmospheric mixing classes; however our data are not sufficient to repeat the analysis for our location.

For Figure 13, the maxima were filtered. For every $x\left(t\right)$, we compute AM9(t) := AM[x(t-4),...,x(t),...,x(t+4)], that is, the 9-hour mean around x(t) leaving out x(t); then, considering that the Poisson SD equals $\sqrt{AM}$, set a 2 sigma threshold $th=AM+2\sqrt{AM}$; and finally, if the hour of the day is between 1 and 8 in the morning and x(t)>th, retaining x(t). One can recognize in Figure 13 that the maxima vary aperiodic considerably, with some weeks with low and some weeks with high maxima.

5.2. Autoregressive Modelling

5.2.1. Simple ARMA

Table 5 shows the AR and AM coefficients of the time series estimated with the ARMA module of Past software. The model has been optimized by the by minimizing the AIC (Akaike information criterion). Figure 14 show the original, modelled and residual series for Vienna (residual = original - model). Perhaps 5 AR components would not really be necessary, as already an ARMA(1,1) model has AIC=2.589e4, that is, the additional autoregressive coefficients improve the model only little. Introducing a differencing step, ARIMA(p,1,1), does not lead to improvement. On the other hand, the difference series x(t)-x(t+24), which largely removes the diurnal cycle, leads to ARMA(1,1) with lower AIC=2.816e4. The model performance may be compared between the series by computing the mean square of the relative residual, defined here as $AM{[res/((original+model)/2)-1]}^2$. We see that the models for the indoor series perform similarly, the one for the outdoor series a bit worse. In the absence of comparable results from the literature we cannot interpret the result in relation to Rn series recorded elsewhere.

5.2.2. Partial Autocorrelation

The shape of the partial autocorrelation function (PACF) informs about the order of the autoregressive model as the number of significant lags. Figure 15 shows the ACF and PACF of the Viennese series. The number of significant PACF lags, 4 to 5, coincides with above findings (Table 5), corroborated also with Gretl software. First order differencing or introducing an ARIMA order d>0 does not improve the model but suggests the presence of one AM component, as also found above (Table 5). For the series: Berlin-out, Berlin-in1 and Berlin-in2, one estimates 4, 7 to 8 and 3 AR components, respectively, slightly different from the findings above. The 7 to 8 AR components suggested for Berlin-in1 from the PACF are also not confirmed by the ARMA analyses of Gretl which suggests 1 to 2 AR components.

5.3. Fractal Analysis

Several authors have applied fractal methods to explore the dynamic of Rn time series. While estimation of the Hurst exponent has led to consistent results, in most cases, other statistics have been less studied and often led to questionable results. For example, [84] applied fractal methods to Rn emissions from a mud volcano in Italy, demonstrating methodology but without conclusive results.

5.3.1. Hurst Exponent

The R/S graph from which the Hurst exponent is derived by log-log linear fitting, is shown in Figure 16 as example for the Vienna series.

In a previous study [57], for the same apartment in Vienna, H=0.92 was found, using an adapted R/S method. The same study reports 0.74-0.92 for different rooms of a house in Salzburg. For outdoor Rn, [85] found H=0.81. A similar list is given in Donner [49]. In a number of studies mainly on indoor and soil Rn series, similar H values were found, indicating persistence and long memory, as could be expected by weather episodes that last typically some days and which are the drivers of aperiodic Rn variability. As three more examples, [86] found H=0.70 – 0.85 for indoor Rn series, [71], H=0.88 and [87], up to 0.9. For Rn exhalation from uranium tailings, [88] found H=0.80. In a study aimed to search for seismic signals in soil Rn time series, [89] found H=0.7-1.

However, there are also studies cited in [85] that found H<0.5, for example [90], opposed to our own study and the ones quoted above. This is implausible on physical grounds and could be caused by erroneous software. Older versions of CDA, such as the freely available v.1.1 from 1995, appear to yield systematically low H values. While a more recent version (v.2.2) is documented in [74], its high cost has limited testing and verification of potential improvements. CDA v.1.1. can be run under DosBox on Window 7 but certain functions are very slow. Table 6 shows Hurst coefficients of the four series estimated with different software and algorithms. Finally, the function ’hurstexp’ in the R package ’pracma’ [91] was used which gives different outputs based on different methods. For both, a large minimum box size was specified. Two are given in Table 6. In spite of a consistent tendency, the differences show that exact estimation seems not to be trivial. Some issues of estimation of the Hurst exponent are discussed in [92].

5.3.2. Multifractal $f\left(\alpha \right)$ Spectra

The $f\left(\alpha \right)$ spectra of the four Rn series are shown in Figure 17. The parabolic shape is as expected; slight deviation from strictly parabolic may point to deviation from multi-fractality. The aperture of the parabolas is different, with the outdoor series showing the highest spread between ${\alpha }_{min}$ and ${\alpha }_{max}$. This points to particularly pronounced "anomalies" (maxima) while the indoor series are comparatively "smoother". This is plausible from a physics perspective because the indoor space serves as a smoother of the dynamic. So far we are not aware of $f\left(\alpha \right)$ spectra of Rn series reported in the literature, therefore we cannot compare our results with the ones of other authors.

5.3.3. Graph Dimension

The box dimensions D sub B of the four time series are listed in Table 6 together with the derived Hurst exponents. The values coincide more or less; in particular for the Berlin outdoor series we recognize a discrepancy whose source is not known. The graphs from which the DB are determined by regression are shown in Figure 18. The values were calculated with the sliding box method, where rectangles of sizes corresponding to the scale are moved over the 2-dim graph space and the rectangles counted which contain points of the graph, that is the number of non-empty boxes which cover the space.

5.3.4. Garger Exponent

The graph from which the exponents g are determined by linear regression is shown in Figure 19. The used range is $1\le T\le 80$ h. The resulting exponents (Table 7) are 0.66 - 0.78 for the indoor and 0.79 for the outdoor series. The exponent of the Viennese series is similar to the one found in [57] for the same room years ago. The series fluctuate more strongly than expected for turbulent flow (g=0.33), but less than white noise (g=1). Outdoor Rn seems to be noisier than indoor Rn, which is physically plausible.

5.3.5. Attractor Embedding Dimension and Lyapunov Exponent

According to Eckmann [93] the highest minimal E that can be reasonably estimated equals approximately 2 lg N, with N the length of the time series. For the Viennese, Berlin outdoor and Berlin indoor-1 series, this implies $E\le 7$, for the shorter series Berlin indoor 2, $E\le 6$.

With CDA 1.1 software, the minimal embedding dimension E can be estimated with the plateau of the correlation dimension versus E or with the false nearest neighbour (FNN) statistic. Delay times $\tau$ were chosen according to the ACF, i.e., the time after which the ACF is close to zero (Table 8). Additionally a standard $\tau$=16 h was used; the results are not very different. They are summarized in Table 8 for $\tau$ according to the ACF. Note that CDA allows setting $\tau$ only in powers of 2. For comparison, in [57], E = 4 to 6 were found. CDA has a module for the Lyapunov exponent, but this does not work under CDA in DosBox. Screen shots of CDA outputs are shown in Figure 20 for Berlin indoor 1 and Vienna.

With R, packages ‘tseriesChaos’ [61] and ‘nonlinearTseries’ [62], we found the following results. Some of the used functions require extensive parameter setting, which was partially done by trial and error.

In the package ’tseriesChaos’, E was estimated with the function ’false.nearest’, the maximal Lyapunov exponent L with ’$lya{p}_k$’ and ’lyap’. For the delay $\tau$, the values according to the ACF close to zero and the value where the first minimum of the ACF occurs were used. Results are not too different. In the package ’nonlinearTseries’, the ’estimateEmbeddingDim’ was used, based on the Cao method [94]; but no plausible or interpretable results could be found, mostly because E was estimated too high, that is $E>7$. For the Lyapunov exponent, the ‘maxLyapunov’ function was used. In all cases, model parameter such as minimal search vicinity and Theiler window (to exclude spurious neighbours) were varied until some plausible result was achieved. (This does not appear very satisfying, but seems to be common practice when analyzing real-world time series.) The function does not give uncertainties or confidence intervals of L, nor p-values to test against the null hypothesis L=0; but from visual impression of the S versus t graphs (S – the separation of phase space points, t – the time along the trajectories from some starting point; graphs not shown here) the values given in Table 9 do not appear significant.

Evidently there is a discrepancy between L values computed with different programs. The low maximal Lyapunov exponents indicate low degree of chaotic behaviour of the series. However, estimation is difficult and results should not be taken too literally, in particular in the absence of much literature results to compare.

We also tried the R package ’DChaos’ [95] but this did not appear successful, yielding consistently negative L; however, these seem to by typical for periodic series. Perhaps the algorithm identified this effect. Further, the software GChaos [96] seems to be user friendly and offers good functionality, but lacks documentation. It applies three algorithms for the Lyaponov exponent, with two methods yielding values around 0.005, one method, values ten times higher.

Planicic et al. [90] reported Lyapunov exponents of hourly sampled outdoor Rn data about ten times higher than the ones of ‘tseriesChaos’, but similar to ones found with ‘nonlinearTseries’. Similarly, Glushkov et al. [97] found maximal L = 0.018 for outdoor Rn. Bossew et al. [57] found values between 0.08 and 0.14 for indoor series.

As a conclusion of this section, stable estimation of the Lyapunov exponent of the time series of this study has not been satisfactory and further investigation is required, including diving deeper into the mathematical and computational details of the algorithms.

6. Conclusions

Although hourly measurement values for low Rn concentrations recorded with the RadonEye monitor have high uncertainty and are often below the theoretical detection limit, it seems that the Rn dynamics can be captured reasonably well. For outdoor environmental applications, the use of the RadonEye should be confined to measurements near the surface at inland sites, where radon concentrations are more likely to exceed the detection limit. It is important to note that this instrument is generally unsuitable for absolute radon measurement in most outdoor environmental conditions, if high accuracy and precision are required.

A 24-hour periodicity is observed in all time series, least significantly in the indoor series recorded in Berlin during the winter period. In this case, we hypothesize that the Rn dynamics are more strongly influenced by those of geogenic Rn carried upwards by the stack effect, rather than by outdoor Rn, as has been found in summer.

The measurements reported here can reproduce findings about Rn concentration cycles given in the literature based on experiments using professional and more sophisticated equipment. However, further work with "cheap" and "simple" equipment, concerning (a) better understanding of the statistical behavior of the RadonEye monitor and (b) more sophisticated experimental setup may prove its potential for atmospheric and geophysical research. The main objective of this work was only to explore the essential feasibility of such simple techniques.

Parallel indoor and outdoor measurements in the same building showed a delayed response of indoor Rn concentration of about 6 hours. This value is a mean over a period and may be variable during the period. This should be elucidated in further studies. Furthermore, such delay value is characteristic for a building and will certainly vary between buildings and rooms within buildings - also a matter which deserves further research.

Interpretation of Rn values reported by the RadonEye is less trivial than it may appear. This is due to the low count rates resulting from low Rn concentrations, and from one feature of its algorithm, namely rounding to integers. The relationship between recorded and estimated true Rn concentration is not trivial. The problem has been approached by Bayesian inversion. An important issue on our wish list for the manufacturer is information about the algorithm implemented in the RadonEye.

Apart from traditional statistical analyses which did not yield particularly surprising or unexpected results, we performed certain fractal analyses of the time series. Except for the Hurst exponent, hardly any literature data are available for comparison. The Hurst exponents were found to be 0.64 for the outdoor and 0.76 - 0.84 for the indoor Rn series, which is consistent with literature. The values point to "persistence" and "long memory" structures, meaning that there are episodes of low or high Rn levels which last for some time. Most probably, this is physically caused by the same behavior of its physical controls, that is, meteorological episodes. Embedding dimensions of the attractors are typically 4 to 6, which may be understood as an estimate of the numbers of relevant control variables. Estimation of Lyapunov exponents turned out to be difficult and yielded partly contradictory results. One reason may be the low degree of chaoticity of the Rn series. Testing against the null hypothesis "no deterministic chaos present" was beyond the scope of this study; it would probably imply extensive simulation with shuffled data.

For further investigation, we propose three main objectives:

• Examine the relationship between recorded dynamics and their controlling physical environmental factors, particularly meteorological variability, with results intended for future publication.
• Refine the statistical methodology through ARFIMA (Autoregressive Fractionally Integrated Moving Average) analysis, leveraging the identified long-term memory structures; enhanced seasonal decomposition; application of detrended fluctuation analysis to reveal fractal structures; and development of techniques to identify chaotic patterns.
• Continue long-term measurements at the established locations to gain deeper insights into the temporal dynamics of radon concentrations and potential patterns that may emerge from extended datasets.