How reliable are process-based ²²²radon emission maps? Results from an atmospheric ²²²radon inversion in Europe

The radioactive noble gas radon (²²²Rn) is a suitable tracer for atmospheric transport and mixing processes that can be used to evaluate and calibrate atmospheric transport models or to estimate greenhouse gas (GHG) emissions using the so-called radon tracer method (RTM). However, these applications require reliable estimates of the ²²²Rn fluxes from the soil. This study evaluates two process-based ²²²Rn flux maps in central Europe in 2021 using the flux results from a 1-year ²²²Rn inversion. The maps are based on different soil moisture reanalysis products (GLDAS-Noah and ERA5-Land), which are used to describe the diffusive ²²²Rn transport in the soil. The ²²²Rn inversion was conducted using the CarboScope-Regional inversion system and observational data from 17 atmospheric sites in central Europe in 2021. We observe that, in particular, the ERA5-Land-based ²²²Rn flux map underestimates the data-driven fluxes from the inversion. Our inversion yields ca. 20 % (GLDAS-Noah) to almost 100 % (ERA5-Land) larger ²²²Rn fluxes than the respective process-based prior fluxes within a domain covering Germany. Also, the temporal variability seems to be underestimated by the process-based flux maps. Using a flat (uniform) prior inversion, we found a significant anti-correlation of −0.6 (and −0.8) between the posterior ²²²Rn flux and the GLDAS-Noah (and ERA5-Land) soil moisture time series, indicating that soil moisture is an important driver for the temporal variability in the ²²²Rn fluxes. To investigate the robustness of our flux estimates, we run the inversion with three different transport models (STILT, FLEXPART, NAME). The respective annual mean flux results agree within ca. 10 %.

1 Introduction

Inverse modelling (e.g. Newsam and Enting, 1988) is a well-established method for constraining surface fluxes of greenhouse gases (GHGs) by minimizing the mismatch between observed and simulated atmospheric dry air mole fractions of the respective GHG. However, limitations in atmospheric transport models, such as inadequate description of vertical mixing within the planetary boundary layer (PBL), can lead to systematic biases in such top-down flux estimates (e.g. Schuh et al., 2019, 2022; Munassar et al., 2023). Therefore, careful quantification of transport model uncertainties is essential for a reliable flux estimation. As this is far from straightforward, potential systematic biases in atmospheric transport models are often assessed by using an ensemble of different transport models in the inversion framework to investigate their impact on the flux estimates (e.g. Geels et al., 2007; Peylin et al., 2013; Monteil et al., 2020; Schuh et al., 2022) or by employing computationally expensive ensemble methods (e.g. Steiner et al., 2024).

A more direct way of quantifying transport model uncertainties is to compare modelled and measured atmospheric activity concentrations of the radioactive noble gas radon (²²²Rn), which is the first gaseous component in the decay chain of uranium that escapes from the soil into the atmosphere (Karstens et al., 2015). As its lifetime (3.8 d) is comparable to the ventilation timescale of the PBL, atmospheric ²²²Rn observations over land contain suitable information on vertical mixing (Jacob and Prather, 1990) and can be used to validate (and possibly even improve) the performance of atmospheric transport models in this respect (e.g. Jacob and Prather, 1990; Chevillard et al., 2002; Gupta et al., 2004; Galmarini, 2006; Zhang et al., 2008; Arnold et al., 2010; Zhang et al., 2021). However, this requires accurate ²²²Rn flux maps that are suitable for modelling atmospheric ²²²Rn activity concentrations.

There are various global and regional ²²²Rn flux maps available, which differ in the methods and the complexity used to describe the ²²²Rn exhalation from the soil (e.g. Rasch et al., 2000; Conen and Robertson, 2002; Zhuo et al., 2008; Szegvary et al., 2009; Griffiths et al., 2010; López-Coto et al., 2013; Karstens et al., 2015; Karstens and Levin, 2024). In an extensive study, Karstens et al. (2015) developed two process-based ²²²Rn flux maps for Europe and compared them with existing ²²²Rn flux maps from the literature. Their study revealed large spatiotemporal differences between the different ²²²Rn flux maps, which can be on the order of the ²²²Rn fluxes themselves, and illustrates the substantial uncertainty associated with ²²²Rn flux maps.

In principle, continuous ²²²Rn flux measurements can be used to validate (and calibrate) the ²²²Rn flux maps (Griffiths et al., 2010; Manohar et al., 2013; Karstens et al., 2015). However, such measurements are sparse and typically only representative for very local spatial scales and often contradict large-scale flux estimates due to the high degree of the soil parameter inhomogeneities. In a technical report on their ²²²Rn flux model, Karstens and Levin (2024) compared the modelled ²²²Rn flux with continuous ²²²Rn flux measurements from one site during 6 months and showed that it is difficult to validate the ²²²Rn flux map of an entire continent with a sparse set of flux measurements. Therefore, long-term, continuous, and reliable ²²²Rn flux measurements with a high spatial representativeness across different parts of Europe are needed for further validation studies. For this reason, Karstens and Levin (2024) recommended performing a ²²²Rn inversion to evaluate the ²²²Rn flux maps, which is a more representative approach, as the atmosphere integrates fluxes from larger areas. In our study, we implemented their suggestion and investigated whether we could evaluate the quality of ²²²Rn flux maps using inverse modelling.

However, there is a conflict in performing a ²²²Rn inversion, since it intrinsically assumes that atmospheric transport is well known and without systematic biases (Schuh et al., 2019). As this is not the case (see above), there is a risk that the inversion will adjust the ²²²Rn fluxes to compensate for unknown biases in the transport model (especially if the transport model uncertainties are not correctly described in the inversion framework). For example, if the vertical mixing in the transport model were too strong, the model would underestimate the ²²²Rn activity concentrations, even if the ²²²Rn fluxes are correct. As a consequence, the inversion would falsely increase the ²²²Rn fluxes. Thus, such incorrectly adjusted ²²²Rn fluxes are useless for modelling ²²²Rn concentrations and validating other transport models.

Therefore, we use an ensemble of transport models and an ensemble of a priori flux maps in our inversion system to carefully quantify the impact of potential biases in the transport model and in the prior fluxes on the ²²²Rn inversion results. These transport models differ, for example, in the parameterization of turbulent motion and convection and in the underlying meteorological data (see Sect. 2.3). Thus, by analysing the influence of the transport models on the ²²²Rn flux estimates, we can assess the robustness of our inversion results. Furthermore, we only use afternoon observations (or nighttime observations for mountain sites), when the atmosphere is typically well mixed and the models are expected to perform best (Gerbig et al., 2008). By doing so, we aim to obtain reliable ²²²Rn flux estimates that can be used to evaluate process-based ²²²Rn flux maps.

In our study, we evaluate the two process-based ²²²Rn flux maps from Karstens and Levin (2022a, b), which have been updated and further refined within the framework of the traceRadon project (Röttger et al., 2021). The maps are based on two different soil moisture products used to describe the ²²²Rn transport through the soil (see Sect. 2.1). We investigate if we can use our ²²²Rn inversion results to assess which of the two ²²²Rn flux maps is best suited for a domain in central Europe well covered by ²²²Rn observations and if we can derive realistic estimates of their uncertainties (Sect. 3.2). Furthermore, we want to investigate whether we can learn something about the ²²²Rn exhalation process from the inversion results, e.g. what information about soil moisture variability is contained in the ²²²Rn observations (Sect. 3.3). In this context, we will also try to improve the process-based ²²²Rn flux maps by using different soil moisture data (Sect. 3.4).

2 Methods

Figure 1 provides an overview sketch of the ²²²Rn inversion process. The individual compartments are described in the following subsections.

Figure 1 Sketch of the inversion process.

2.1 Process-based ²²²Rn flux maps

In this study, we investigate different process-based ²²²Rn flux maps, all of which are based on the ²²²Rn flux model developed by Karstens et al. (2015) for an “infinitely deep unsaturated homogeneous soil”. We will mainly focus on the ²²²Rn flux maps from Karstens and Levin (2022a, b), which are based on GLDAS-Noah and ERA5-Land soil moisture and porosity data to describe the ²²²Rn transport through soil. However, we also want to analyse a third, alternative ²²²Rn flux map that we have prepared for this study using high-resolution soil moisture and porosity data from the German Weather Service (DWD; see Table 1). In the following, we briefly describe this ²²²Rn flux model, but the reader is referred to Karstens et al. (2015) and Karstens and Levin (2024) for more comprehensive details.

Table 1 Overview of the ²²²Rn flux maps used in this study.

The main assumption of the ²²²Rn model is that ²²²Rn exhalation from the soil occurs mainly by diffusive transport (Nazaroff, 1992). This assumption is justified because the ²²²Rn activity concentration in the air in the soil is several orders of magnitude higher than in the ambient air above the soil (Čeliković et al., 2022). Assuming steady-state conditions and a ²²²Rn source Q that is constant with depth, the following equation can be derived for the ²²²Rn flux j at the soil surface (z=0):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}j(z=\mathrm{0})& =-Q\sqrt{{\displaystyle \frac{D_{\mathrm{e}}}{\mathit{\lambda }}}}\tanh \left(z_G\sqrt{{\displaystyle \frac{\mathit{\lambda }}{D_{\mathrm{e}}}}}\right)\\ & =-Q\stackrel{\mathrm{‾}}{z}\tanh \left({\displaystyle \frac{z_G}{\stackrel{\mathrm{‾}}{z}}}\right).\end{array}\end{array}$$

Hence, the flux j (in Bq m⁻² s⁻¹) depends on the ²²²Rn source Q (in Bq m⁻³ s⁻¹), the effective diffusion coefficient D_e (in m² s⁻¹), the ²²²Rn decay constant λ (in s⁻¹), and the water table depth z_G (in m). The ²²²Rn relaxation depth $\stackrel{\mathrm{‾}}{z}$ (in m) is given by $\stackrel{\mathrm{‾}}{z}=\sqrt{D_{\mathrm{e}}/\mathit{\lambda }}$. The ²²²Rn source Q can be described as a product of the concentration of radium (²²⁶Ra, i.e. the precursor of ²²²Rn) in the soil, the ²²²Rn decay constant, the dry bulk density of the soil, and the emanation coefficient, which describes the probability that the ²²²Rn atoms can escape from the soil grains, in which they were formed, into the soil air. According to the parameterization from Zhuo et al. (2008), this emanation coefficient depends on the soil texture, soil moisture, soil porosity, and soil temperature. The ²²⁶Ra concentration was calculated using a map of the uranium (²³⁸U) content in the soil from the European Atlas for Natural Radiation (EANR; Cinelli et al., 2019) and assuming secular equilibrium between ²³⁸U and its daughter ²²⁶Ra (see Karstens and Levin, 2024).

The effective diffusion coefficient D_e is described with the parameterization from Millington and Quirk (1960), which depends on soil porosity and soil moisture. Karstens and Levin (2024) used the average soil moisture and porosity data from the upper 40 cm of the soil to calculate the diffusion coefficient (and the emanation coefficient). Hence, they assumed that the ²²²Rn released from the soil surface originates mainly from the top 40 cm of the soil. In Sect. 3.4, we revisit this assumption. The ²²²Rn flux model also takes into account a temperature dependence of the diffusion coefficient according to the parameterization from Schery and Wasiolek (1998).

Furthermore, the ²²²Rn flux from the soil is reduced in regions with shallow water table depths as the soil water hinders the ²²²Rn exhalation. It is described by the factor $\tanh \left(\frac{z_G}{\stackrel{\mathrm{‾}}{z}}\right)$ in Eq. (1). For water table depths z_G that are large compared to the ²²²Rn relaxation depth $\stackrel{\mathrm{‾}}{z}$ (i.e. $z_G\gg \stackrel{\mathrm{‾}}{z}$), this factor becomes 1. All three ²²²Rn flux maps used in this study were calculated using Eq. (1) but with different soil moisture data products as listed in Table 1. The two ²²²Rn flux maps from Karstens and Levin (2022a, b), based on GLDAS-Noah and ERA5-Land soil moisture, respectively, cover a slightly smaller area compared to our model domain. Therefore, we have spatially extended these ²²²Rn flux maps by filling the extended (land) areas with the respective daily mean land fluxes of the whole domain (see Fig. 3a, b). In the following, we refer to these two ²²²Rn flux maps as “GLDAS” and “ERA5”, respectively.

For the third ²²²Rn flux map, we make use of the high-resolution soil moisture data from the AgrarMeteorologische Berechnung der Aktuellen Verdunstung (AMBAV; Herbst et al., 2021) model, which is a water balance model operated by the German Weather Service (DWD) for agricultural purposes that provides daily, 1 km² resolution soil moisture data for Germany. Information on regional soils from detailed geological maps (Hartmann et al., 2024) is thereby incorporated into the model. The AMBAV model calculates soil moisture data for crop and grass-covered soil types. Whilst the GLDAS-Noah and ERA5-Land models assume that the soil porosity is constant over the entire soil column, AMBAV uses vertically resolved soil porosity data (from Hartmann et al., 2024). We combined the AMBAV soil moisture product with soil moisture data for forests simulated with the forest hydrological model LWF-Brook90 (Hammel and Kennel, 2001) to obtain a high-resolution soil moisture map for the dominant land cover types in Germany, which we refer to as “DWD” soil moisture hereinafter. The basis for the classification of the pixels into arable land, grassland, and forest is the land cover data from Corine Landcover (2018), the classification of the forest pixels into the main tree species (beech, oak, spruce, and pine) was carried out using the tree species map by Blickensdörfer et al. (2022). To be consistent with Karstens and Levin (2024), we initially use the average DWD soil moisture (and porosity) of the top 40 cm of the soil. However, in Sect. 3.4, we further investigate how the vertical averaging of the soil moisture impacts the ²²²Rn fluxes. Since the DWD soil moisture is only available for Germany, we fill the area outside Germany (and the few gaps over German cities) with the GLDAS-Noah soil moisture data so that we can calculate a complete ²²²Rn flux map for Europe. To be consistent, we also apply the porosity data used by the GLDAS-Noah model for areas outside Germany. We call this third ²²²Rn flux map “DWD/GLDAS”. All three ²²²Rn flux maps are re-gridded to a horizontal resolution of 0.05° × 0.05°.

2.2 Atmospheric ²²²Rn concentration observations

We use hourly ²²²Rn activity concentration observations from 17 European sites in 2021. The ²²²Rn measurements have been performed with several detector types, which are based on different measurement principles and assumptions. Figure 2 and Table 2 give an overview of the different observation sites and the detectors used. The main characteristics of the used detector types are compiled in Schmithüsen et al. (2017) and Grossi et al. (2020). For example, the radon detectors developed by the Australian Nuclear Science and Technology Organisation (ANSTO) are based on the so-called dual-flow-loop two-filter approach, which provides a direct measure of the ²²²Rn activity concentration (Griffiths et al., 2016). In this method, the sampled air passes through an initial filter that removes all ambient aerosols and ²²²Rn progenies. The filtered air is then directed into a large delay volume (typically > 1000 L) in which new ²²²Rn progenies (²¹⁸Po and ²¹⁴Po) are formed. A second flow loop within the delay volume frequently circulates the sample of air through a second filter to ensure that all ²²²Rn progenies are collected on this filter (e.g. ²¹⁸Po has a half-life of only 3 min) and that their α decay can be counted. Finally, the ²²²Rn activity concentration is calculated from the α decay of the ²²²Rn progenies and the flow rate. Due to the large delay volume, ANSTO detectors have a relatively slow response time of about 45 min. Approximately 40 % of the signal is observed 1 h after the radon pulse is delivered (Griffiths et al., 2016). This delay necessitates time response correction. To recover the true atmospheric signal, the ANSTO data from the sites in UK are corrected using a deconvolution procedure, following the methodology described in Griffiths et al. (2016) and further validated and standardized in Kikaj et al. (2025). This approach effectively recovers rapid variations in ²²²Rn activity concentrations, which are typically observed during periods of changing air mass origin, such as transitions between terrestrial and oceanic fetches or during the early morning and late evening hours when the planetary boundary layer (PBL) height changes rapidly. For the non-UK ANSTO data, we applied a first-order deconvolution by shifting the ²²²Rn data by 1 h. Since only hourly averaged data from the afternoon period are used in this study, when atmospheric mixing is strongest and ²²²Rn concentrations are relatively stable, such a first-order deconvolution might be acceptable for our purposes. Furthermore, ANSTO detectors are affected by a gradually increasing background signal due to the accumulation of long-lived ²¹⁰Po in the measurement chamber. This background is routinely monitored every 2–3 months and extracted from the signal using standardized correction procedures as part of the regular data processing protocol, ensuring the accuracy of reported ²²²Rn activity concentrations (Röttger et al., 2025; Kikaj et al., 2025).

Figure 2 Map of the European model domain. The Germany domain is indicated by the grey rectangle. The locations of the sites with ²²²Rn observations used in the inversion are marked by coloured dots. The colours indicate the different radon detectors described in the text.

Table 2 Overview of the 17 European sites providing atmospheric ²²²Rn measurements.

T – tower sites; M – mountain sites; S – coastal sites; U – urban sites. ^a High-altitude sites for which a correction height was used in the model to account for the steep terrain.
^b Sites for which the FLEXPART/NAME particle release height is > 20 m different than that for STILT.

In contrast, the Heidelberg Radon Monitor (HRM; Levin et al., 2002; Gachkivskyi and Levin, 2022) detectors installed at the German sites measure the ²²²Rn activity concentrations indirectly. In the atmosphere, the ²²²Rn progenies get attached to aerosols. The HRM detectors collect these atmospheric aerosols on a filter and measure the α decay of the ²²²Rn progenies. In order to determine the ²²²Rn activity concentration from those measurements, assumptions must be made about the radioactive disequilibrium between ²²²Rn and its progenies (Schmithüsen et al., 2017). This disequilibrium depends on the height above ground (Jacobi and André, 1963). Intercomparison studies between ANSTO and HRM detectors revealed that radioactive equilibrium between ²²²Rn and its daughter products is reached between ca. 50–100 m above ground (Schmithüsen et al., 2017; Grossi et al., 2020). We therefore applied equilibrium correction factors for observation sites with an air intake height below 90 m above ground. Furthermore, wet deposition of atmospheric aerosols and aerosol loss in long sampling lines can lead to artefacts in the HRM-based ²²²Rn activity concentrations (Xia et al., 2010; Grossi et al., 2016; Levin et al., 2017). To account for this, we (1) flagged the HRM data during situations with high air humidity > 98 % (> 95 % for the mountain sites HPB, SSL, and TOH) as suggested by Gachkivskyi et al. (2025) and (2) applied an aerosol loss correction, as described in Levin et al. (2017), for sites with sampling tubing lengths > 15 m. Finally, we averaged the half-hourly HRM ²²²Rn data to obtain hourly ²²²Rn observations. The corrected and moisture-selected HRM ²²²Rn observations are compiled in Fischer et al. (2024).

The radon detector installed at the Mace Head (MHD) observation site is another type of one-filter radon monitor, again based on an indirect method of determining ²²²Rn activity concentration. It was developed at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE) in France (Biraud, 2000) and uses a moving filter band system to collect and measure the ²²²Rn progenies ²¹⁸Po and ²¹⁴Po. Since Schmithüsen et al. (2017) found a ²¹⁴Po $/$ ²²²Rn disequilibrium factor of 1.0 for this coastal site, we did not apply a disequilibrium correction to the MHD measurements. However, we again applied a humidity flag to account for the wet deposition of atmospheric aerosols. Note that the MHD ²²²Rn observations have a temporal resolution of only 2 h.

Finally, we converted all ²²²Rn observations to the HRM detector scale in order to obtain a harmonized data set. Note that the ANSTO detector scale would result in approximately 11 % larger ²²²Rn activity concentrations compared to the HRM detector scale (Schmithüsen et al., 2017). The LSCE measurements from MHD must be divided by 0.95 to convert them to the HRM scale (Schmithüsen et al., 2017). We assume an uncertainty of 0.5 Bq m⁻³ for the hourly ²²²Rn observations used in the inversion to account for instrumental uncertainties, uncertainties associated with the required calibrations and corrections (Grossi et al., 2020), and uncertainties in the background contributions from outside Europe, which we subtract from the ²²²Rn observations (see Sect. 2.3). As mentioned above, we only use observations from PBL sites during well-mixed situations in the afternoon (11:00–16:00 UTC), when the transport model is expected to show the best performance. In contrast, nighttime (23:00–04:00 UTC) observations are used for mountain sites when the impact of local (thermally induced) wind systems is expected to be negligible.

Note that there was a leak in the ²²²Rn inlet line at the Cabauw (CBW) site in 2021, which may have contaminated the ²²²Rn measurements at this site. However, as we found no obvious anomalies when comparing the observation–model differences at CBW with those at the nearby Lutjewad (LUT) site, we decided to use the CBW data in our study. Nevertheless, we investigated the impact of the CBW observations on the inversion results, which turned out to be small (see Appendix C).

2.3 Transport models

For our main inversion runs, we use the Stochastic Time-Inverted Lagrangian Transport model (STILT; Lin et al., 2003). In Lagrangian transport models, an ensemble of numerical particles is released from each site every time step (here every hour of the measurements), and their back-trajectories are computed by confronting the particles with the underlying meteorological fields and a stochastic representation of turbulent motions. From the distribution of the back-trajectories, so-called footprints are deduced, which describe the sensitivity of the observation site to surface fluxes at the individual pixels in the catchment area of the site. The ²²²Rn activity concentrations C(x_s,t_s) at the observation site x_s and time t_s can then be calculated by convolving the ²²²Rn flux F(x_j,t_i) of each grid cell x_j (j=1, …, N, with N being the number of grid cells) and time step t_i (i=1, …, T, with T being the maximal duration of the back-trajectories, i.e. 240 h in our case) with the footprint f and taking into account the radioactive decay of ²²²Rn with a lifetime of τ_Rn=5.52 d:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}C\left(x_s,t_s\right)=& \sum _j^N\sum _i^{T}f(x_s,t_s|x_j,t_i)\cdot F\left(x_j,t_i\right)\\ & \cdot \exp \left(-{\displaystyle \frac{t_i}{{\mathit{\tau }}_{\mathrm{Rn}}}}\right).\end{array}\end{array}$$

STILT is driven with meteorological fields from the Integrated Forecasting System of the European Centre for Medium-Range Weather Forecasts (ECMWF-IFS), extracted at a spatial resolution of 0.25° × 0.25° (lat × long) and a temporal resolution of 3 h. Every hour, 100 particles are released from the observation sites, and their back-trajectories are calculated for 10 d or until they leave the STILT model domain shown in Fig. 2. In STILT, the horizontal resolution of the footprint matrix is dynamically reduced as the spatial distribution of the STILT particles increases. This prevents the under-sampling of surface fluxes at times and regions where the STILT particles are distributed over extensive areas and have large gaps between each other (Gerbig et al., 2003). Therefore, this enables the usage of a small number of only 100 particles in STILT and thus reduces computational costs while ensuring proper statistics. The final STILT footprints are mapped on a grid with a horizontal resolution of 0.25° × 0.25°.

Table 3 Overview of the different transport models.

In order to assess the influence of the transport model on the inversion results, we perform two additional inversion runs using footprints from two alternative Lagrangian transport models (see Table 3): the FLEXible PARTicle dispersion model (FLEXPART; Stohl et al., 2005; Pisso et al., 2019) and the Numerical Atmospheric-dispersion Modelling Environment (NAME; Jones et al., 2007). Similarly to STILT, FLEXPART is also driven by meteorological analysis and forecasts from the operational ECMWF-IFS high-resolution (HRES) runs, available hourly at 0.1° × 0.1° resolution for the European domain and 3-hourly at 0.5° × 0.5° for the rest of the globe. With FLEXPART, 20 000 particles are released every hour, and the back-trajectories are calculated for 30 d or when leaving a domain encompassing parts of North America, the North Atlantic, and Europe. The resulting footprints are mapped on a grid with a horizontal resolution of 0.234° × 0.352° (lat × long). NAME utilizes meteorological data from the UK Met Office Unified Model (UM; Cullen, 1993). As with FLEXPART, 20 000 particles are released every hour, and the back-trajectories are calculated for 30 d. The NAME footprints also have a horizontal resolution of 0.234° × 0.352°. The FLEXPART and NAME footprints have initially been calculated for another project and another model domain. Therefore, the back-trajectories were calculated for 30 d, to ensure that most of the particles have left this larger model domain. However, to resolve ²²²Rn fluxes in Europe, a 10 d simulation, as is done with STILT, is sufficient.

Moreover, the FLEXPART and NAME footprints that we are using as alternative footprints for this study are already post-processed, meaning that they are (1) corrected for the ²²²Rn decay and (2) aggregated over the duration of the back-trajectories (i.e. over 30 d or until the particles have left the model domain). This means that we have no temporal information about the particle distributions any more. However, since the ²²²Rn flux maps have a daily resolution, we need to know for which ²²²Rn flux time period the time-aggregated footprints are mainly sensitive. To investigate this, we make use of the time-resolved STILT footprints and perform two types of ²²²Rn simulations for three different sites (continental, coastal, and mountain sites; see Fig. A1). In the first simulation, we modelled ²²²Rn concentrations by mapping the time-resolved STILT footprints with the daily ²²²Rn fluxes (as described in Eq. 2). For the second type of simulation, we firstly aggregated the STILT footprints over the duration of the back-trajectories (i.e. over 10 d) and corrected them for the ²²²Rn decay to mimic the post-processed FLEXPART and NAME footprint products. We then mapped these time-aggregated STILT footprints with ²²²Rn fluxes, which were averaged over different time intervals ranging from 1 to 10 d (see Fig. A1). By comparing both types of simulations, we can investigate for which ²²²Rn flux time interval the time-aggregated footprints are most sensitive. We found that the ²²²Rn concentration differences between the first and the second simulations show a low bias and a small standard deviation when the time-aggregated footprints are mapped with ²²²Rn fluxes that are averaged over 3 d, indicating that the time-aggregated footprints are mainly representative for the average ²²²Rn fluxes of the first 3 d before particle release. Therefore, we average the ²²²Rn fluxes over 3 d when mapping them with the FLEXPART and NAME footprints. However, we want to emphasize that we use the time-resolved STILT footprints to perform our main analyses in this study and that we apply the time-aggregated FLEXPART and NAME footprints only to investigate the impact of different transport models on the inversion results, by which means we assess the robustness of our results.

The back-trajectories of the numerical particles are only calculated whilst the particles remain within the European model domain. For particles that leave the model domain, lateral boundary conditions are needed. For this purpose, we have constructed a simple global ²²²Rn flux map, assuming a constant ²²²Rn flux of 1 atom cm⁻² s⁻¹ (21 mBq m⁻² s⁻¹) over land surfaces south of 60° N, a halved ²²²Rn flux of 0.5 atoms cm⁻² s⁻¹ (10.5 mBq m⁻² s⁻¹) over land in the higher latitudes > 60° N, and a vanishing ²²²Rn flux in permafrost regions. The much smaller ²²²Rn fluxes from the ocean were neglected. We use the Eulerian global atmospheric Tracer Model (TM3; Heimann and Körner, 2003) and this global ²²²Rn flux map to simulate hourly ²²²Rn activity concentrations for each European grid cell. A ²²²Rn background concentration is then determined for each site and hour by averaging the modelled ²²²Rn activity concentrations in the respective grid cells of the endpoints of the STILT back-trajectories after taking radioactive decay into account. Due to the relatively short ²²²Rn atmospheric lifetime, these background ²²²Rn activity concentrations are typically quite low (the median of the mean background contribution across all sites is < 15 %). We subtract this modelled ²²²Rn background from the ²²²Rn observations and use the resulting ²²²Rn excess concentrations to constrain the ²²²Rn fluxes in Europe. Note that we use the same ²²²Rn excess concentrations when we apply the time-aggregated FLEXPART and NAME footprints in the inversion.

To investigate the impact of the simulated background on the ²²²Rn flux results, we perform two additional inversion runs using alternative ²²²Rn backgrounds. For the first run, we replace the TM3-simulated ²²²Rn concentration field by the ²²²Rn concentration field provided by the Copernicus Atmosphere Monitoring Service (CAMS) global reanalysis (EAC4; Inness et al., 2019) and again use the endpoints of the STILT trajectories to calculate an alternative (decay-corrected) ²²²Rn background. For the second run, we fully neglect the Rn background; i.e. we assume that it is zero.

2.4 Inversion setup

In this study, we use the CarboScope-Regional (CSR) inversion system described in Rödenbeck et al. (2003, 2009) to constrain the process-based ²²²Rn fluxes with the atmospheric ²²²Rn observations. In the following, we provide a brief overview of the CSR inversion system, with a specific focus on the aspects that are particularly relevant to the ²²²Rn inversion. For more technical details about the inversion algorithm and how the iterative solution is found, we refer the reader to Rödenbeck (2005).

In the CSR system, the Bayesian cost function is implemented using a “linear statistical flux model” (Rödenbeck et al., 2009):

$$\begin{array}{}\text{(3)}& \mathit{x}={\mathit{x}}_{\mathrm{fix}}+\mathbf{F}\mathit{p}.\end{array}$$

The flux field x is written as the sum of a fixed flux component x_fix, which is the mean of the prior flux field x_prior, and the product between a matrix F and a vector p with adjustable parameters. By construction, the a priori parameters in p have zero mean and unit variance and are uncorrelated. The columns of F represent spatiotemporal flux patterns, which are scaled by the elements of p. The prior error covariance matrix B is then given by B=FF^T, and the cost function J(p) is

$$\begin{array}{}\text{(4)}& J\left(\mathit{p}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{p}}^{\mathrm{T}}\mathit{p}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{c}}_{\mathrm{obs}}-\mathbf{H}\mathit{x}\right)}^{\mathrm{T}}{\mathbf{R}}^{-\mathrm{1}}({\mathit{c}}_{\mathrm{obs}}-\mathbf{H}\mathit{x}).\end{array}$$

The second term in Eq. (4) shows the data constraint, and the first term describes the prior constraint.

The quadratic Bayesian cost function is minimized by applying a conjugate gradient algorithm that allows large state vectors. The cost function includes a model–data mismatch vector that contains the hourly differences between the observed (c_obs) and modelled (c_mod=Hx, with H being the transport matrix) ²²²Rn activity concentrations from all 17 sites. The model–data mismatch vector is weighted with a covariance matrix R containing the uncertainties of the ²²²Rn observations and the transport model. As already mentioned, we assume an uncertainty of 0.5 Bq m⁻³ for the hourly ²²²Rn observations. The uncertainty of the transport model is chosen dependent on the type of observation site (Rödenbeck 2005). For example, continental tower sites can typically be better represented in the model than sites such as Heidelberg, which is located in a narrow river valley with complex local circulation. Depending on the site, we therefore assume a transport model uncertainty of ca. 0.5 to 1.5 Bq m⁻³ (see Table 2). The total model–data mismatch uncertainty is obtained by adding the observation and transport model uncertainties quadratically. To account for temporal correlations between model–data mismatch within 1 week, we applied the so-called data density weighting proposed by Rödenbeck (2005). This inflates the uncertainty of the model–data mismatch by the square root of the number of observations within a week. The typical timescale for synoptic events is 1 week, in which the ²²²Rn observations are expected to be temporally correlated. As mentioned above, we only use afternoon observations (or nighttime observations for mountain sites) when the atmosphere is typically well mixed. In addition, we applied a 2σ filtering to these data (as described in Rödenbeck et al., 2018) to exclude the observations with the largest model–data mismatch, as these are considered to be inadequately represented by the transport model.

The Bayesian approach adds a priori information to stabilize the solution. We use the process-based ²²²Rn flux maps described in Sect. 2.1 and a map with spatially and temporally constant fluxes over the European continent (“flat prior”) as a priori estimates. Due to the large differences between the ²²²Rn flux maps, we assume an a priori uncertainty of 100 % for the European ²²²Rn fluxes and for a timescale similar to the temporal correlation length. In the standard setting of our inversion system, we assume that the a priori flux errors are spatially correlated over a length scale of about 400 km, and we choose a temporal correlation length of 3.5 d (“Filt52T” in CarboScope notation). In Sect. 3.2.2, we investigate how changes to these settings affect the inversion results. The a priori uncertainties are described by the a priori covariance matrix B, which determines the a priori constraint. The ratio between a priori and data constraint determines how strongly the solution is regularized by the a priori information. The inversion system minimizes the model–data mismatch by scaling the ²²²Rn fluxes, taking into account the model–data mismatch and the a priori uncertainties. Overall, we use the CSR system to determine daily ²²²Rn fluxes for the whole year 2021 at a spatial resolution of 0.25°.

2.5 Overview of the different inversion runs

To investigate how robust the ²²²Rn flux estimates are, we perform several inversion runs using different prior flux estimates, varied parameter settings in the inversion setup, and different transport models (see an overview in Table 4).

Table 4 Overview of the different inversion runs.

As soil moisture controls temporal changes in the diffusion coefficient within the soil and is therefore expected to be the main driver of the temporal variability in the ²²²Rn flux, we want to investigate what information about soil moisture is contained in the ²²²Rn observations. For this purpose, we constructed a “flat” prior with spatially and temporally constant ²²²Rn fluxes over land. Such a flat-prior inversion does not use any a priori information on soil moisture variability. This means that the a posteriori flux variability is only caused by the signals in the ²²²Rn observations (possibly with spurious contributions from variations in the transport model error). To investigate the extent to which variations in the ²²²Rn flux are caused by changes in soil moisture, we will calculate the temporal correlation between the daily a posteriori flux of a domain covering Germany and the daily GLDAS-Noah and ERA5-Land soil moisture average of the same domain, respectively (see Sect. 3.3).

3 Results

3.1 Comparison of the process-based ²²²Rn flux maps

Figure 3 shows a comparison between the three process-based ²²²Rn flux maps. The GLDAS and ERA5 ²²²Rn flux maps show similar ²²²Rn hotspot regions, e.g. on the Iberian Peninsula and in Italy, which can be explained by the high uranium activity concentrations there (see Fig. 3a, b). Note that these ²²²Rn flux maps differ only in the soil moisture and porosity data but not in the underlying uranium map. Compared to the ERA5-Land soil moisture, the GLDAS-Noah soil moisture leads to larger annual mean ²²²Rn fluxes in central Europe but to smaller ²²²Rn fluxes in Scandinavia (Fig. 3c). Overall, the annual mean flux differences between both maps can be as large as the fluxes themselves. The DWD/GLDAS ²²²Rn flux map has a higher spatial resolution for Germany (Fig. 3e). The DWD soil moisture and porosity lead to lower annual mean ²²²Rn fluxes in the southern part of Germany and to slightly higher ²²²Rn fluxes in northern Germany than the ²²²Rn fluxes based on GLDAS-Noah soil moisture and porosity data (Fig. 3f).

Figure 3 Annual mean ²²²Rn fluxes in 2021 for the European model domain (a–b) and for the Germany domain (d–e) based on GLDAS-Noah (a, d), ERA5-Land (b), and DWD/GLDAS (e) soil moisture data. The DWD/GLDAS ²²²Rn flux map (e) is based on DWD soil moisture and porosity data for Germany and on GLDAS-Noah soil moisture and porosity data for outside Germany. Panels (c) and (f) show the annual mean differences between the GLDAS-Noah- and ERA5-Land-based ²²²Rn fluxes and the GLDAS-Noah- and DWD/GLDAS-based ²²²Rn fluxes, respectively.

3.2 Top-down evaluation of the GLDAS and ERA5 ²²²Rn flux maps

3.2.1 Comparison between a posteriori and process-based ²²²Rn fluxes

We start by evaluating the existing ²²²Rn flux maps from Karstens and Levin (2022a, b), which are based on the GLDAS-Noah and ERA5-Land soil moisture data. On an annual average, the inversion increases the process-based a priori ²²²Rn fluxes in central Europe and decreases them over the British Isles (see Fig. 4a, b). This corresponds to the negative model–data mismatches caused by the a priori fluxes at most of the observation sites in central Europe and the positive model–data mismatches at the Irish site Mace Head (MHD) and the British site Weybourne (WAO; see Fig. B1). The negative flux adjustments over the British Isles occur throughout the whole year (see Fig. D1). They could be due to boundary effects, as the British Isles are most affected by potential biases in the ²²²Rn background due to the prevailing westerly winds in Europe. The ERA5 a priori fluxes lead to more negative model–data mismatches at most of the continental sites than the GLDAS a priori fluxes (Fig. B1). Consequently, the ERA5 fluxes in central Europe are increased much more by the inversion than the GLDAS ²²²Rn fluxes, bringing the posterior fluxes closer to each other than the priors (Fig. 4c). While in the case of the GLDAS inversion there are also some negative flux adjustments in France during summer, the ERA5 inversion shows positive flux adjustments in France and Germany throughout the whole year (Fig. D1).

Figure 4 Results of the CSR-STILT inversion runs with GLDAS (red), ERA5 (blue), and flat (grey) a priori ²²²Rn fluxes. Panels (a) and (b) show the annual mean innovation (a posteriori minus a priori ²²²Rn flux differences) for the inversion runs based on the prior fluxes from GLDAS and ERA5, respectively. Panels (c) and (d) show the time series for the full year of the a priori (dotted line) and the a posteriori (solid line) ²²²Rn fluxes averaged over the entire European domain (c) and over the Germany domain (d), respectively, for 2021. The Germany domain is depicted by the magenta rectangle in the maps in panels (a) and (b).

The posterior flux estimate of the flat-prior inversion shows some substantial differences of up to roughly 50 % for individual months compared to the posterior estimates based on the GLDAS and ERA5 prior fluxes (compare the grey curve with the red and blue curves in Fig. 4c). In particular, the European flux estimates only show a substantial seasonal cycle when the process-based prior fluxes are used, which in turn also exhibit a pronounced seasonal cycle. In contrast, using the flat prior results in much weaker seasonal variability in the posterior flux. This demonstrates that the continental-scale ²²²Rn flux estimates and their seasonal variability are strongly influenced by the prior information. This can be explained by the sparse ²²²Rn data coverage in Europe, which is insufficient to reliably constrain the seasonal variability in the European ²²²Rn fluxes. In the following, we will therefore focus our analysis on a central European domain around Germany, which is covered well by the ²²²Rn observations.

In this Germany domain, the flat-prior inversion yields very similar flux estimates to the inversions based on the GLDAS and ERA5 a priori fluxes (Fig. 4d). The annual means of the three different posterior fluxes agree within 10 %, which is rather small considering the two process-based a priori fluxes in the Germany domain differ by more than 60 % on average in 2021. Also, the temporal variability in the posterior fluxes is comparable; the normalized standard deviations of the three different posterior fluxes agree within 10 %. Furthermore, the different inversion runs produce similar seasonal variations, which are comparable with the seasonal variation in the process-based prior fluxes. While the soil-moisture-driven GLDAS and ERA5 prior fluxes are 29 % and 37 % higher in the summer half-year of 2021 than in the winter half-year, the corresponding posterior fluxes show 36 % and 27 % higher fluxes in summer than in winter, respectively. The flat-prior inversion leads to a very similar result: the flux estimates in the summer half-year are 30 % higher than in the winter half-year. Overall, this illustrates that the flux estimates in the Germany domain are indeed well constrained by the observations and less affected by the prior information. Therefore, in the following, we will use the flat-prior inversion results, which are not based on additional information about temporal and spatial flux variability (and are thus based on the least a priori information), to evaluate the performance of the process-based ²²²Rn flux maps.

In the Germany domain, the annual mean a posteriori ²²²Rn flux (in the flat-prior inversion) is about 20 % and almost 100 % larger than the GLDAS and ERA5 process-based ²²²Rn fluxes, respectively, indicating that, in particular, the ERA5 ²²²Rn fluxes might be too low in central Europe. In addition, the a posteriori flux shows a higher temporal variability than the bottom-up fluxes. The normalized standard deviation of the daily a posteriori ²²²Rn flux is about 10 % (ERA5) and 35 % (GLDAS) higher than the normalized standard deviation of the process-based fluxes in the Germany domain. In Sect. 3.4, we investigate if we can increase the process-based ²²²Rn fluxes and their temporal variability by using higher-resolution soil moisture and porosity data for Germany. Overall, the inversion reduces the annual mean and the standard deviation of the model–data mismatch at almost all sites (Fig. B1).

3.2.2 Sensitivity of inversion results to model configuration

After having shown that the different a priori fluxes (GLDAS, ERA5, flat) lead to only small changes in the a posteriori fluxes in the Germany domain, we want to investigate how robust these inversion results are. For this, we conduct several inversion runs using different settings, to assess sensitivity (see Table 4).

Figure 5 (a) Daily a posteriori ²²²Rn fluxes in the Germany domain for different inversion runs performed with CSR-STILT. Panel (b) shows the relative difference between the respective a posteriori fluxes and the flat-prior inversion run (grey curve in panel a) with the standard parameter settings described in Sect. 2.4. The different inversion runs are described in Table 4.

The various inversion runs performed with CSR-STILT lead to very similar a posteriori ²²²Rn fluxes in Germany (see Fig. 5). These runs differ in terms of the prior information used, the inversion's uncertainty parameters, the scale of the ²²²Rn observations, and the simulated ²²²Rn background. Most of these sensitivity tests yield a posteriori fluxes with an annual mean difference well below 10 % compared to the a posteriori flux of the flat-prior inversion run based on the standard setting (described in the first row of Table 4). This shows that we indeed get robust results for the Germany domain, which is well covered by ²²²Rn observation sites. There is one exception to this result, however: the inversion run “flat_ANSTO”, which yields an a posteriori flux that has an almost constant offset of 13 % compared to the respective inversion run with the standard setting (magenta curve in Fig. 5). For this run we have applied a different ²²²Rn observation scale (ANSTO scale instead of HRM scale; see Sect. 2.2; note that the simulated background ²²²Rn concentration is the same as for the inversions with the HRM scale). Since this ANSTO scale leads to 11 % higher ²²²Rn activity concentrations than the HRM scale, such an offset in the ²²²Rn fluxes is to be expected. This illustrates the urgent need for well-calibrated and SI-traceable ²²²Rn observation data sets.

Next, we analyse the results of the CSR-FLEX and CSR-NAME inversions, for which FLEXPART and NAME footprints are used instead of the STILT footprints. When the GLDAS prior ²²²Rn fluxes are mapped with the NAME footprints, the resulting ²²²Rn concentrations strongly overestimate the observations from the mountain site Schauinsland (SSL) by ca. 2.7 Bq m⁻³ on an annual average (see Fig. B1). Such large overestimations are not observed at other sites and in the STILT and FLEXPART simulations. This could be due to the fact that the NAME model, unlike STILT and FLEXPART, did not use an elevated release height for the SSL site (see Table 2) and therefore has difficulty in correctly representing the steep terrain there. Therefore, we excluded the observations from SSL in the case of the CSR-NAME inversion to avoid unrealistic flux adjustments in southwestern Germany.

Figure 6 Comparison between the CSR-STILT (grey), CSR-FLEX (green), and CSR-NAME (pink) inversion results for the Germany domain. All inversions were performed with the flat a priori flux (grey dotted line in panel a). In the case of the CSR-NAME inversion, the observations from SSL were not used (“w/o SSL”) due to unexpectedly large prior model–data mismatches at that site (see text). Panel (a) shows the daily ²²²Rn fluxes, and panel (b) shows the relative difference between the different a posteriori fluxes and the flat-prior inversion run performed with CSR-STILT (solid grey curve in panel a).

For most of the continental sites, FLEXPART and NAME show higher surface influences and thus larger prior ²²²Rn concentrations and smaller model–data mismatches than STILT (Fig. B1). Averaged over all sites in the Germany domain, the FLEXPART simulations lead to 46 % and 18 % smaller model–data mismatches than STILT (for GLDAS and ERA5 ²²²Rn flux, respectively), and the NAME simulations lead to 47 % and 32 % smaller model–data mismatches than STILT (average of NAME model–data mismatch without SSL). Consequently, the CSR-FLEX and CSR-NAME inversions lead to smaller ²²²Rn flux adjustments in central Europe than STILT. On an annual average, the posterior ²²²Rn flux based on the flat prior is 5 % and 12 % smaller for the CSR-FLEX and CSR-NAME inversions, respectively, than for the CSR-STILT inversions in the Germany domain (see Fig. 6). However, there seems to be a slight seasonal cycle in the difference. In the winter half-year, the FLEXPART- and NAME-based a posteriori ²²²Rn fluxes are within 3 % of the STILT-based fluxes, whereas, during summer, the FLEXPART- and NAME-based a posteriori ²²²Rn fluxes are 12 % and 23 % lower than the STILT-based flux, respectively. This seasonal cycle in the difference could indicate seasonal variations in the model performances, which should be investigated further. As a consequence, the FLEXPART and NAME runs exhibit a smaller seasonal cycle in the ²²²Rn flux than the STILT run.

Overall, these results show that changing the transport model leads to some deviations in the a posteriori fluxes, but these are comparable to the deviations caused by varying the inversion parameter settings, as shown in Fig. 5. Moreover, these transport-model-based deviations are not larger than the annual flux differences caused by the choice of the ²²²Rn observation scale.

3.3 What soil moisture information is in the ²²²Rn observations?

Figure 7 shows the temporal correlation between the daily a posteriori ²²²Rn flux and soil moisture data, both averaged over the Germany domain. For both soil moisture products, we obtain quite strong anti-correlations between about −0.6 for GLDAS-Noah and even −0.8 for ERA5-Land. This anti-correlation means that we get high ²²²Rn fluxes when soil moisture is low, which makes sense, as soil pores are less filled with water in dry conditions and the diffusion coefficient of ²²²Rn is higher in air than in water. In addition, the stronger anti-correlation in the case of ERA5-Land might indicate that this reanalysis product describes the temporal variability in the soil moisture better than GLDAS-Noah. We also interpret the existence of meaningful correlations as a confirmation that the inversion is indeed picking up real flux variations (even though spurious correlations due to the relationships of soil moisture and atmospheric mixing cannot be excluded either).

Figure 7 (a) Flat prior (in grey, dotted) and associated a posteriori flux (in grey, solid) for the Germany domain, together with the ERA5 and GLDAS a posteriori fluxes for comparison. (b) ERA5-Land (blue) and GLDAS-Noah (red) soil moisture within the Germany domain for the upper 0–10 cm (solid), 0–40 cm (dashed), and 0–100 cm (dotted) of the soil column. (c) Correlation between the posterior flux curve of the flat-prior inversion shown in panel (a) and the different soil moisture curves shown in panel (b), calculated for different time lags of the soil moisture curves.

3.4 Evaluating soil depth sensitivity in ²²²Rn flux

An open question is why the process-based flux maps driven by soil moisture fields show less variations than the flux estimates from the inversion. Karstens et al. (2015) estimated that the soil parameters of the top 100 cm of soil are most important for ²²²Rn flux at the surface, and Karstens and Levin (2024) used the average soil moisture of the top 40 cm of soil in their ²²²Rn flux model. To investigate systematically over which soil depth interval the soil moisture data should be averaged, we calculated the average soil moisture of the top 100 cm, top 40 cm, and top 10 cm of soil for both reanalysis products GLDAS-Noah and ERA5-Land, and we calculated time-lagged correlations with the a posteriori flux of the flat-prior inversion (Fig. 7c). The average soil moisture of the top 10 cm shows the largest temporal variability and increases rapidly when it rains. The increase in the average soil moisture of the top 40 cm and top 100 cm is slightly delayed after a rain event (cf. Figs. 7 and G1).

If we average the soil moisture data of the top 40 cm of soil, as done in Karstens and Levin (2024), the maximum anti-correlation between the soil moisture time series and the posterior ²²²Rn flux in the flat-prior inversion is reached when the soil moisture curve is shifted by 1 d, which means that the soil moisture lags the ²²²Rn flux by 1 d. This time lag even increases to 2–3 d if the soil moisture of the top 100 cm of the soil is averaged. Therefore, the ²²²Rn flux responds faster to rain or drought events than the average soil moisture of the top 40 or 100 cm of soil, suggesting that the average soil moisture responds in a delayed way. However, if only the soil moisture of the top 10 cm is averaged, the time lag between soil moisture and ²²²Rn flux disappears. This could indicate that the variability in the ²²²Rn flux is mainly caused by the soil moisture variability in the top 10 cm of the soil. This leads to the following question: does the ²²²Rn flux model produce larger temporal variations if driven by average soil moisture data from the top 10 cm of soil only?

Figure 8 Soil moisture (a), air-filled pore space (b), and process-based ²²²Rn fluxes (c) based on the ERA5-Land (blue), GLDAS-Noah (red), and DWD/GLDAS (orange) soil moisture and porosity data averaged over the top 10 cm (solid) and the top 40 cm (dashed) of the soil column. The a posteriori flux (grey) based on the flat prior is shown in panel (c) for comparison.

Indeed, the average soil moisture in the top 10 cm of the soil shows a higher temporal variability than the average soil moisture in the top 40 cm, leading to an increased temporal variability in the ²²²Rn flux (see Fig. 8). Since the average soil moisture of the top 10 cm is higher during rain events and lower during dry periods than the average soil moisture of the top 40 cm, the ²²²Rn fluxes based on the 10 cm average soil moisture tend to be higher during dry periods and lower during rain events compared to the ²²²Rn fluxes based on the 40 cm average soil moisture (see Fig. G1). The daily ²²²Rn fluxes can differ by up to roughly 30 % (ERA5) and 60 % (GLDAS) if the 10 cm soil moisture is used instead of the 40 cm soil moisture. The normalized standard deviation of the ²²²Rn fluxes based on the top 10 cm soil moisture data is 22 % (ERA5) and 24 % (GLDAS) higher compared to the respective fluxes based on the top 40 cm soil moisture data. Consequently, the bottom-up ²²²Rn fluxes based on the top 10 cm soil moisture data show a similar temporal variability to the flux estimate of the flat-prior inversion (the normalized standard deviations of the bottom-up fluxes and the inversion result agree within about 10 %). Nevertheless, using the top 10 cm soil moisture data instead of 40 cm has, for both reanalysis products, only a very minor impact on the annual mean ²²²Rn fluxes (in fact, the soil moisture data from the top 10 cm of soil lead to ca. 2 %–3 % lower annual mean ²²²Rn fluxes). It should be noted that both reanalysis products use the same porosity values for all soil layers, so they assume that porosity does not change with soil depth.

Next, we evaluated if we could further improve the process-based ²²²Rn fluxes by using a high-resolution soil moisture data product from the DWD model for Germany, which is based on vertically resolved porosity data (see Sect. 2.1). The average DWD soil moisture in the top 10 cm lies between the ERA5-Land and GLDAS-Noah soil moisture and is on average even 0.02 m³ m⁻³ higher than in the top 40 cm. However, the porosity in the top 10 cm is also higher than in the top 40 cm, which could be explained by a layer of humus below the soil surface. Thus, the air-filled pore space, which results from the porosity minus the soil moisture, is on an annual average very similar in the top 10 cm and in the top 40 cm of the soil. Moreover, the air-filled pore space is not larger for the DWD data product than for GLDAS-Noah. Therefore, even with the high-resolution soil moisture and the vertically resolved porosity data from the DWD model, the ²²²Rn fluxes are still roughly 25 % smaller compared to the a posteriori flux in the Germany domain. However, the bottom-up ²²²Rn fluxes based on the DWD data also show a higher (15 %) temporal variability for the 0–10 cm soil moisture average than for the 0–40 cm average. Hence, the normalized standard deviation of the DWD-based ²²²Rn flux based on the top 10 cm soil moisture data is again similar to the normalized standard deviation of the a posteriori ²²²Rn flux.

4 Discussion

4.1 Does the ²²²Rn data coverage allow robust flux estimates in central Europe?

So far, inverse modelling has only rarely been used to estimate ²²²Rn fluxes. We are aware of a study by Hirao et al. (2010), who estimated the Asian ²²²Rn fluxes with a Bayesian inversion using atmospheric ²²²Rn observations from seven sites in eastern Asia and an Eulerian transport model. Their results indicate a higher ²²²Rn flux in eastern Asia than suggested by an a priori estimate from a ²²²Rn exhalation model. However, to our knowledge, a ²²²Rn inversion has not yet been performed over Europe. Therefore, we firstly investigated whether the current coverage of atmospheric ²²²Rn observations is sufficient to estimate robust ²²²Rn fluxes over Europe. We found that the annual mean and the temporal variability in the ²²²Rn flux estimates for the whole European model domain strongly depend on the prior information used. This can be explained by the poor coverage of available ²²²Rn observations in large parts of Europe. Currently, the Integrated Carbon Observation System (ICOS) atmospheric station network is planning to release hourly resolved ²²²Rn activity concentration measurements from several tower sites on a regular basis. Provided the number of ICOS sites carrying out ²²²Rn measurements increases in the future, as was suggested previously (ICOS RI, 2020), this will improve the availability of ²²²Rn data in Europe.

Our inversion results are much less affected by the prior information if we restrict the domain to a region around Germany, a region well covered by ²²²Rn observations. For this region, the differences between the annual mean ²²²Rn flux estimates of various inversion configurations are much smaller (on the order of 10 %) compared to the annual mean differences between the process-based ²²²Rn fluxes (> 60 %). This indicates that the current data coverage in central Europe enables robust inversion results that are mainly determined by the observational data and that are relatively independent of the choice of the a priori fluxes. This is a prerequisite for reliably evaluating process-based ²²²Rn fluxes with an atmospheric transport inversion.

The fact that our ²²²Rn observation data set cannot constrain the high-flux regions in Spain, Portugal, or Italy could affect the inversion results for our regional target domain in central Europe, if fluxes from these regions are transported to Germany by advective winds. To assess this effect, we performed another inversion run using additional ²²²Rn observations from the mountain site Puy de Dôme (PUY), which is located in the southern part of France, where the process-based ²²²Rn flux maps indicate elevated fluxes, too (see Fig. 3). The impact of the PUY site on the annual mean ²²²Rn flux in the Germany domain is less than 3 % (see Appendix E). As the high-flux regions on the Iberian Peninsula and in Italy are further away from our target region than the PUY site, the effects of dilution and radioactive decay should be even larger. Therefore, we expect these high-flux regions to have only a minor influence on the ²²²Rn flux estimates for the Germany domain.

Finally, beyond the importance of good data coverage in the region of interest, differences in the scale of the ²²²Rn observations directly translate into the inversion results and determine the extent to which absolute ²²²Rn fluxes can be estimated with a ²²²Rn inversion (see Fig. 9). This illustrates the need for further intercomparison projects between the ²²²Rn instruments and a well-calibrated and SI-traceable ²²²Rn data set. A proposed standardized protocol for the harmonization of ²²²Rn measurements has recently been published by Kikaj et al. (2025) thanks to the traceRadon project.

Figure 9 Annual mean ²²²Rn fluxes in the Germany domain for 2021. Shown are the process-based GLDAS, ERA5, and DWD/GLDAS prior ²²²Rn fluxes and the a posteriori fluxes based on the respective prior fluxes. The hatched range indicates the increase in the a posteriori fluxes when the ANSTO ²²²Rn observation scale is used instead of the HRM scale.

4.2 What is the performance of process-based European ²²²Rn flux maps?

Our inversion leads to ca. 20 % (GLDAS) and almost 100 % (ERA5) larger annual mean ²²²Rn fluxes than the process-based ²²²Rn fluxes in the Germany domain, and the temporal variability is also higher in the posterior flux. Estimating correct absolute ²²²Rn fluxes with a ²²²Rn inversion is challenging due to the aforementioned differences in the scale of the ²²²Rn observations and due to general uncertainties of atmospheric transport inversions, such as the representation of vertical mixing or lateral boundary conditions. In fact, we found that the FLEXPART and NAME transport models lead to up to 12 % lower annual mean ²²²Rn fluxes in comparison with STILT. In contrast, using the ANSTO scale instead of the HRM scale would lead to even higher ²²²Rn flux estimates and thus further increase the bias compared to the process-based fluxes. The effect of the lateral boundary conditions on the annual mean ²²²Rn fluxes in central Europe is on the order of a few percent: using alternative boundary conditions based on CAMS instead of TM3 results in 6 % higher annual mean fluxes in the Germany domain, and using zero boundary conditions would lead to even 12 % higher fluxes (see Fig. F1). From that, we conclude that especially the ERA5 ²²²Rn fluxes might be underestimated in central Europe and that the bias compared to the posterior flux is unlikely to be fully explained by deficits in the inversion. This is consistent with other studies, which found that the ERA5 soil moisture might be too high in central Europe (Li et al., 2020; Karstens and Levin, 2024).

The significant (anti-)correlation between the posterior flux of the flat-prior inversion and the soil moisture data indicates that the temporal variability in the inversion signals is reliable and not only caused by inversion noise. Compared to GLDAS-Noah, the ERA5-Land soil moisture leads to a larger anti-correlation, indicating that it may better describe the temporal variability in the soil moisture in central Europe. This finding seems to be corroborated by a direct comparison with soil moisture measurements from several sites in Germany (see Fig. H1). Overall, the differences between the a priori and the a posteriori ²²²Rn fluxes might give a rough estimate for the uncertainty of the process-based ²²²Rn fluxes (see Fig. 9).

4.3 What can we learn from the inversion results to improve the ²²²Rn flux maps?

The comparison between the temporal variability in the posterior ²²²Rn flux estimate of the flat-prior inversion and the soil moisture variability shows that soil moisture is an important driver for the temporal variability in the ²²²Rn fluxes and that the ²²²Rn flux variability is mainly driven by soil moisture in the top soil layer. The latter could be explained by the following points: (1) ²²²Rn gas produced in deeper soil layers has to travel a longer distance in the soil before it reaches the atmosphere and may be more affected by radioactive decay than ²²²Rn gas produced in upper soil layers. The soil moisture in the upper soil layers would therefore have a greater effect on the ²²²Rn flux at the surface than soil moisture in the deeper soil layers. (2) Depending on the type of soil, precipitation may take days to reach the deeper layers of the soil (cf. Fig. G1). As a result, air from the soil layers closer to the surface will outgas first, and the outgassing of air from deeper layers will be delayed. Soil moisture averaging over several layers may dilute the expected correlation between soil moisture and ²²²Rn flux at the surface. (3) Similarly, if the upper layer is saturated after smaller amounts of precipitation, but there are still many air-filled pores underneath, the top saturated layers reduce the air-filled pore space and therefore diffusion potential. ²²²Rn from deeper soil layers will hardly reach the surface.

The temporal variability in the process-based ²²²Rn fluxes increases for all three soil moisture products, GLDAS, ERA5, and DWD, when ²²²Rn fluxes are calculated using the soil moisture data from the top 10 cm instead of the top 40 cm of the soil column, resulting in ²²²Rn fluxes with similar temporal variability to our a posteriori flux. This finding could also open up the possibility of creating ²²²Rn flux maps using satellite-based soil moisture retrievals, which are only sensitive to the top few centimetres of soil. Even though such satellite data can have large gaps, e.g. due to snow cover, frozen ground, or strong vegetation cover, such an application for improved ²²²Rn flux maps might be possible.

Absolute ²²²Rn fluxes are strongly dependent on the air-filled pore space, which determines the strength of the diffusion process. The comparison between GLDAS-Noah and ERA5-Land illustrates that an annual mean difference of 0.06 m³ m⁻³ in the air-filled pore space already leads to a difference of more than 60 % in the annual mean ²²²Rn fluxes in the Germany domain. However, using the soil moisture and porosity data from the top 10 cm instead of the top 40 cm of soil has almost no effect on the annual mean ²²²Rn fluxes. This is because, in the case of GLDAS-Noah and ERA5-Land, the annual mean soil moisture is similar for the top 10 cm and the top 40 cm, and the porosity is vertically constant. The DWD soil moisture is higher in the top 10 cm than in the top 40 cm. However, this is compensated by the (vertically resolved) porosity, which is also higher in the top 10 cm than in the top 40 cm. Hence, the annual mean of the air-filled pore space and of the ²²²Rn flux are similar for the top 10 cm and top 40 cm soil data. If a vertically averaged porosity is assumed for the DWD data (i.e. if the average porosity of the top 40 cm is also used to calculate the air-filled pore space of the top 10 cm of soil), the air-filled pore space of the top 10 cm would be ca. 0.02 m³ m⁻³ smaller than that of the top 40 cm. Consequently, the resulting ²²²Rn fluxes based on the top 10 cm instead of the top 40 cm soil moisture would be on the order of 10 %–20 % lower. Therefore, soil moisture impacts not only the temporal variability in the ²²²Rn fluxes, but, together with porosity, also the absolute ²²²Rn fluxes. Thus, an overestimation of the soil moisture and/or an underestimation of the porosity can easily change the absolute ²²²Rn fluxes. Overall, a reduction in the differences in the absolute values of different soil moisture and porosity products would further constrain process-based ²²²Rn fluxes. In this respect, an expansion of soil moisture measurement networks, such as the International Soil Moisture Network (ISMN), the Cosmic-ray soil moisture monitoring network (COSMOS), or similar, would be desirable. Furthermore, cosmic-ray neutron sensing provides soil moisture data with a spatial representativeness of a few hundred metres (Köhli et al., 2015). This might make it better suited to evaluating (coarse-resolution) soil moisture reanalysis products, which can show significant differences compared to local (point) measurements (cf. Fig. H1).

Of course, there may also be other parameters or processes that are not or only inadequately described in the ²²²Rn flux model and could explain a bias in the absolute fluxes, e.g. an underestimation of the ²²²Rn source due to radium concentrations or the emanation coefficient being too low. In addition, advective fluxes, e.g. induced by a pressure gradient between the atmosphere and the soil, could lead to a ²²²Rn flux contribution. However, their overall contribution might be negligible compared to the diffusive flux because the permeability of the soil is several orders of magnitude smaller than the diffusion coefficient, and the hydrostatic pressure gradients are typically small (López-Coto et al., 2013; Nazaroff, 1992). Nevertheless, changes in atmospheric pressure, e.g. due to weather fronts, can induce short-term variability in the ²²²Rn flux: decreasing atmospheric pressure is expected to suck ²²²Rn-rich air from the soil into the atmosphere, whereas rising atmospheric pressure forces ²²²Rn-poor air from the atmosphere into the soil (Clements and Wilkening, 1974). However, due to the compensating effects of decreasing and rising atmospheric pressure, the overall effect of pressure is larger on instantaneous ²²²Rn flux variability than on absolute (time-averaged) ²²²Rn fluxes (Schery et al., 1984; Ishimori et al., 2013).

Indeed, we found no significant temporal correlation (R² < 0.05) between the a posteriori ²²²Rn flux and the atmospheric pressure (and its first derivative) in the Germany domain and also for a smaller domain in southwestern Germany. This may indicate that the ²²²Rn flux variability is much more influenced by soil moisture than by pressure changes on the timescale accessible to the inversion (i.e. daily resolution) and/or that the pressure effect is superimposed by other meteorological effects (e.g. wind, precipitation). Overall, our inversion results provide information on ²²²Rn exhalation processes occurring on a daily or longer timescale. The processes driving sub-daily ²²²Rn flux variability, which can be on the order of 20 %–50 % (Rábago et al., 2022), may be difficult to investigate using our inversion results.

4.4 Outlook

Overall, further development of the ²²²Rn flux maps is needed. Our inversion results could already give some indications on how to further improve the maps. It would also be interesting to compare our inversion results with other ²²²Rn flux maps that are based on different ²²²Rn flux models or even completely different methods, such as ²²²Rn flux estimation using the terrestrial γ dose rate (Szegvary et al., 2009). While Szegvary et al. (2009) presented a “static” inventory based on specific and punctual measurements, continuous γ dose rate time series, e.g. from the European Radiological Data Exchange Platform (EURDEP), may enable time-resolved γ-dose-rate-based ²²²Rn flux maps. These maps could then be used alongside our inversion results to further investigate temporal variability in ²²²Rn fluxes.

Finally, accurate ²²²Rn flux maps are important not only for evaluating atmospheric transport models or constraining the emissions of other trace gases, e.g. using the radon tracer method (RTM; e.g. Grossi et al., 2018; Levin et al., 2021; Curcoll et al., 2025). They are also relevant to the radiation protection community because ²²²Rn is one of the major sources of natural radiation to which the population is exposed. Nowadays, European countries are developing national plans for achieving their radiation protection goals. These plans need to count on reliable ²²²Rn flux maps with a high spatial and temporal resolution for assessing so-called radon priority areas (RPAs; Röttger et al., 2021). Therefore, the evaluation of ²²²Rn flux maps using a ²²²Rn inversion might be of particular interest to both the climate and the radon protection community.

5 Conclusions

The characteristics of ²²²Rn make it a powerful natural tracer for atmospheric transport that can be used to validate and calibrate atmospheric transport models or to estimate GHG fluxes, e.g. with the radon tracer method (RTM; e.g. Levin et al., 2021). However, all these applications require an accurate estimate of the ²²²Rn flux. So far, ²²²Rn fluxes have mainly been evaluated by local ²²²Rn flux measurements, which have a limited spatial representativeness. In this study, we investigated the potential of a top-down approach to evaluate process-based ²²²Rn fluxes, although we are well aware that such an approach is cyclic and assumes unbiased model transport, which of course may not be the case (otherwise the information from ²²²Rn observations would no longer be needed). We carefully assessed the impact of potential deficits in the transport models by performing several inversion runs with three different transport models. We found that the current coverage with ²²²Rn observations in Europe allows robust and data-driven ²²²Rn flux estimates within a region covering Germany, which also provide indications on how to improve the process-based ²²²Rn flux maps.

We evaluated two process-based ²²²Rn flux maps for Europe, which are based on two different soil moisture reanalysis products (GLDAS-Noah versus ERA5-Land). Our a posteriori annual mean flux of about 14 ± 4 mBq m⁻² s⁻¹ (mean ± SD) is ca. 20 % (GLDAS) and almost 100 % (ERA5) higher than the process-based ²²²Rn fluxes in a domain covering Germany in 2021. Furthermore, both ²²²Rn flux maps tend to underestimate the temporal variability in the ²²²Rn fluxes, although the variability in the ERA5-based ²²²Rn fluxes is in better agreement with the variability in our inversion results.

We found a significant correlation of r = −0.6 and r = −0.8 between the posterior flux of a flat-prior inversion, which is independent of any prior information about spatial and temporal flux variations, and the GLDAS-Noah and ERA5-Land soil moisture, respectively. Moreover, the soil moisture time series lag a few days behind the ²²²Rn flux if the soil moisture is averaged over a depth that is too large (i.e. > 40 cm). In contrast, the time series of the soil moisture and the ²²²Rn flux are time-synchronous if the soil moisture average of the top 10 cm of soil is used. This indicates that the temporal ²²²Rn flux variability is mainly caused by the soil moisture variability in the top 10 cm of soil. Indeed, we were able to increase the temporal variability in the process-based ²²²Rn fluxes by using soil moisture data from the top 10 cm only, resulting in a temporal variability similar to that of the posterior flux estimate.

Finally, realistic uncertainty estimates for the ²²²Rn flux maps are required when ²²²Rn is used in joint inversions for a targeted tracer such as methane (CH₄). Such a dual-tracer inversion directly incorporates the ²²²Rn information on the transport model performance by exploiting the fact that the transport model error (as part of the model–data mismatch error) is correlated between the targeted tracer (e.g. CH₄) and ²²²Rn. In a subsequent study, we will investigate whether this information can help to improve the top-down flux estimates of the targeted tracer (CH₄).

Appendix A: How to average the ²²²Rn fluxes for the time-aggregated FLEXPART and NAME footprints

Figure A1 Modelled ²²²Rn concentration differences between the simulations based on time-aggregated and time-resolved STILT footprints. The time-resolved STILT footprints are mapped with daily ²²²Rn fluxes (based on GLDAS-Noah soil moisture). The time-aggregated STILT footprints are mapped with GLDAS ²²²Rn fluxes, which were averaged over 1 d (grey), 2 d (red), 3 d (blue), 5 d (orange), and 10 d (magenta). This experiment has been performed to deduce the most appropriate time interval over which the daily ²²²Rn fluxes should be averaged when they are mapped with the time-aggregated FLEXPART and NAME footprints. The results are shown for a continental site (KIT200), a coastal site (WAO15), and a mountain site (SSL450). From this study, we found that averaging the ²²²Rn fluxes over 3 d might lead to a good compromise between low bias and low standard deviation between the time-aggregated and time-resolved STILT runs (see Sect. 2.3).

Appendix B: Fits to observations

Figure B1 Fits to the observations for the CSR-STILT, CSR-FLEXPART, and CSR-NAME setups, using the GLDAS and ERA5 a priori fluxes. Shown are annual mean model–data mismatch (a) and its standard deviation (b) for the a priori fluxes (dotted) and the a posteriori fluxes (solid) for each observation site.

Appendix C: Impact of the observations from Cabauw (CBW)

As a leak in the ²²²Rn inlet line at the CBW site may have contaminated the ²²²Rn measurements in 2021 at this site, we investigated the impact of the CBW observations on the inversion results. For this, we performed an additional (flat-prior) inversion run without using the observations from CBW. The resulting posterior ²²²Rn fluxes are a few percent lower in the surroundings of CBW if the CBW observations are not used. However, averaged over the Germany domain, the annual mean flux is only less than 1 % lower.

Figure C1 Impact of the observations from Cabauw (CBW) on the inversion results. (a) Ratio between the ²²²Rn posterior flux of the inversion without CBW data and the standard inversion with CBW data. (b) Time series for the full year of both posterior ²²²Rn fluxes in the Germany domain.

Appendix D: Seasonal posterior–prior ²²²Rn flux differences

Figure D1 Innovation (posterior–prior ²²²Rn flux differences) for CSR-STILT inversion runs using the GLDAS (red) and ERA5 (blue) prior fluxes, averaged for the first (a), second (b), third (c), and fourth (d) quarter of 2021.

Appendix E: Impact of unconstrained high-flux regions on inversion results in the Germany domain

To assess the impact of distant high-flux regions not covered by our observations on the inversion results in the Germany domain, we made a sensitivity inversion run using the additional observations from the mountain site Puy de Dôme (PUY). This site is located in the southern part of France, for which the process-based ²²²Rn flux maps indicate elevated ²²²Rn fluxes (see Fig. 3).

We initially excluded the PUY site in our inversion system due to the following reasons: the PUY site uses an LSCE ²²²Rn instrument, and an intercomparison between the LSCE and the HRM ²²²Rn instruments at that specific mountain site has never been made. Therefore, it is unclear how to scale the PUY measurements to make them comparable to the HRM scale. Furthermore, the PUY site is challenging to represent in the transport model due to its very steep terrain.

Nevertheless, for this PUY inversion run, we assume a scaling factor of 0.7 to make the PUY measurements comparable to the HRM scale. This factor has been found in Schmithüsen et al. (2017) for the French site Gif-sur-Yvette (GIF), which has a similar ²²²Rn intake height above ground to the PUY site. However, while the GIF site is located 167 m a.s.l., the PUY site is located 1465 m a.s.l., which calls into question the suitability of the assumed scaling factor. Therefore, we will only interpret the PUY inversion results in a qualitative way.

Figure E1 shows the posterior flux difference between the PUY inversion, which uses the additional observations from the PUY site, and our standard inversion (without PUY). The largest flux differences occur southwest of the PUY site in the southern part of France and on the Iberian Peninsula, where our base inversion is poorly constrained due to the lack of observations. The southwesterly orientation of the footprint of the PUY site can be explained by the predominant influence of southwesterly winds in that region. Consequently, the PUY observations mainly contain information on the fluxes in southwestern Europe. Indeed, the impact of the PUY observations on the fluxes in the Germany domain are small. When including the PUY observations, the annual mean flux in the Germany domain changes by less than 3 %.

From that, we expect that the potential influence from high-flux regions in southwestern Europe, e.g. in Spain, Portugal, and Italy, on the inversion results in the Germany domain is negligible. Those regions are even more distant than the PUY site, and the effect of dilution and radioactive decay should therefore be even larger.

Figure E1 Posterior ²²²Rn flux difference between an inversion run using and not using the ²²²Rn observations from the PUY site (marked in yellow). Both inversions use the GLDAS prior fluxes.

Appendix F: Impact of ²²²Rn background on inversion results in the Germany domain

Figure F1 Impact of the European ²²²Rn background on the ²²²Rn flux estimates for the Germany domain. Shown are the results of flat-prior inversions using a TM3-simulated background (solid grey curve), a CAMS (EAC4)-based background (orange), and a zero background (magenta).

Appendix G: Difference between ²²²Rn fluxes based on 0–10 and 0–40 cm soil moisture averages

Figure G1 (a) Mean ERA5-Land (blue), GLDAS-Noah (red), and DWD/GLDAS (yellow) soil moisture in the Germany domain, averaged between 0–10 cm (bold line) and 0–40 cm (dashed line), and mean precipitation in Germany (grey histograms). Note that the soil moisture average corresponds to a slightly larger area (Germany domain) than the precipitation average (Germany only, data taken from DWD, 2024). (b) Difference between the soil moisture 0–10 and 0–40 cm averages for the Germany domain. (c) Differences between the ²²²Rn fluxes based on 0–10 and 0–40 cm soil moisture data.

Appendix H: Comparison of soil moisture reanalysis products with observations

Figure H1 Comparison between the GLDAS-Noah (red), ERA5-Land (blue), and DWD (orange) soil moisture with measurements (black; data taken from ICOS RI, 2024) performed at different ICOS ecosystem sites in Germany. The modelled soil moisture shows the average between 0–10 cm, and the measurements show the mean soil moisture from all sensors within the 0–10 cm soil layer. The first three-letter code of the subplot titles indicates the measurement site, and the second three-letter code refers to the site class: gra – grasslands; cro – croplands; dbf – deciduous broadleaf forests; enf – evergreen needleleaf forests.