Bayesian transdimensional inverse reconstruction of the Fukushima Daiichi caesium 137 release

Abstract

The accident at the Fukushima Daiichi nuclear power plant (NPP) yielded massive and rapidly varying atmospheric radionuclide releases. The assessment of these releases and of the corresponding uncertainties can be performed using inverse modelling methods that combine an atmospheric transport model with a set of observations and have proven to be very effective for this type of problem. In the case of the Fukushima Daiichi NPP, a Bayesian inversion is particularly suitable because it allows errors to be modelled rigorously and a large number of observations of different natures to be assimilated at the same time. More specifically, one of the major sources of uncertainty in the source assessment of the Fukushima Daiichi NPP releases stems from the temporal representation of the source. To obtain a well-time-resolved estimate, we implement a sampling algorithm within a Bayesian framework – the reversible-jump Markov chain Monte Carlo – in order to retrieve the distributions of the magnitude of the Fukushima Daiichi NPP caesium 137 (${}^{\mathrm{137}}$Cs) source as well as its temporal discretization. In addition, we develop Bayesian methods that allow us to combine air concentration and deposition measurements as well as to assess the spatio-temporal information of the air concentration observations in the definition of the observation error matrix.

These methods are applied to the reconstruction of the posterior distributions of the magnitude and temporal evolution of the ${}^{\mathrm{137}}$Cs release. They yield a source estimate between 11 and 24 March as well as an assessment of the uncertainties associated with the observations, the model, and the source estimate. The total reconstructed release activity is estimated to be between 10 and 20 PBq, although it increases when the deposition measurements are taken into account. Finally, the variable discretization of the source term yields an almost hourly profile over certain intervals of high temporal variability, signalling identifiable portions of the source term.

1 Introduction

1.1 The Fukushima Daiichi nuclear accident

On 11 March 2011, an earthquake under the Pacific Ocean off the coast of Japan triggered an extremely destructive tsunami that hit the Japanese coastline about 1 h later, killing 18 000 people. These events led to the automatic shutdown of four Japanese nuclear power plants. However, the submersion of emergency power generators prevented the cooling system from functioning, making the shutdown impossible at the Fukushima Daiichi nuclear power plant (NPP). This yielded a massive release of radionuclides, characterized by several episodes of varying intensity, that lasted for several weeks. The estimation of such a release is difficult and subject to significant uncertainties. Reconstructing the temporal evolution of the caesium 137 (${}^{\mathrm{137}}$Cs) source requires one to establish a highly variable release rate that runs over several hundreds of hours.

Since 2011, several approaches have been proposed to assess the Fukushima Daiichi radionuclide release source terms. For instance, methods based on a simple comparison between observations and simulated predictions have been investigated by , , , , , , , and , and ambitious inverse problem methods have been developed in studies such as , , and . In particular, used air concentration and deposition measurements to assess the source term of ${}^{\mathrm{137}}$Cs by rigorously estimating the errors, and evaluated the source term by inverse modelling with ambient gamma dose measurements. Finally, inverse problem methods based on a full Bayesian formalism have also been used. applied several methods, including Markov chain Monte Carlo (MCMC) algorithms, to estimate the time evolution of the ${}^{\mathrm{137}}$Cs release, and they accompanied this with an estimate of the associated uncertainty. refined the ${}^{\mathrm{137}}$Cs source term with an optimization method based on a Bayesian inference from several measurement databases.

1.2 Bayesian inverse modelling and sampling of a radionuclide source

Bayesian inverse problem methods have been shown to be effective in estimating radionuclide sources . Here, we briefly describe this Bayesian framework for source term reconstruction; a more complete description is available in . The crucial elements defining the control variable vector of the source $\mathit{x}$ are (i) $\ln \mathit{q}$, a vector whose components correspond to the logarithm of the release rates $q_i$ (constant releases at a time interval, such as 1 d or 1 h), and (ii) hyperparameters such as entries of the observation error scale matrix $\mathbf{R}$ or the time windows of the release rates (more details in Sect. ).

The posterior probability density function (pdf) of $\mathit{x}$ is given by Bayes' rule:

1$$p\left(\mathit{x}\right|\mathit{y})={\displaystyle \frac{p\left(\mathit{y}\right|\mathit{x}\left)p\right(\mathit{x})}{p\left(\mathit{y}\right)}}\propto p(\mathit{y}\left|\mathit{x}\right)p\left(\mathit{x}\right)$$ where $\mathit{y}$ is the observation vector. The pdf $p\left(\mathit{y}\right|\mathit{x})$ is the likelihood, which quantifies the fit of the source vector $\mathit{x}$ to the data. More precisely, the observations $\mathit{y}$ are compared to a set of modelled concentrations: $\mathbf{H}\mathit{q}$ (i.e. the predictions, constructed using a numerical simulation of the radionuclide transport from the source). $\mathbf{H}$ is the observation operator, the matrix representing the integration in time (i.e. resolvent) of the atmospheric transport model. Therefore, the predictions of the model are considered linear in $\mathit{q}$. The likelihood is often chosen as a distribution parameterized with an observation error scale matrix $\mathbf{R}$. The pdf $p\left(\mathit{x}\right)$ represents the information available on the source before data assimilation.

Once the likelihood and the prior are defined, the posterior distribution can be estimated using a sampling algorithm. Sampling techniques include, in particular, the very popular Markov chain Monte Carlo (MCMC) methods that assess the posterior pdf of the source, or marginal pdfs thereof. This method has been applied by to estimate the Algeciras incident source location, by , and by , who evaluated Xenon-133 releases at Chalk River Laboratories. More recently, the technique has been used to assess the uncertainties in the source reconstruction by , , and .

1.3 The transdimensional analysis and objectives of this study

The problem of finding a proper representation of $\mathit{x}$ (i.e. a discretization of the source term) is the problem of finding the adequate step function for $\mathit{q}$, i.e. the optimal number of steps, the optimal time length of each step, and the corresponding release rates. This representation depends on the case under study: depending on the observations, more or less information is available to define a more or less temporally resolved source term. In other words, depending on the data available to inform one about the source term, the steps of the function supporting $\mathit{q}$ could be small (a lot of available information, allowing for a fine representation) or large (little information available, providing a coarse representation), yielding an irregular temporal discretization of the source term. Note that the choice of the discretization can be seen as a balancing issue due to the bias–variance trade-off principle : low-complexity models are prone to high bias but low variance, whereas high-complexity models are prone to low bias but high variance. The difficulty in the source term discretization is that of selecting representations that are sufficiently rich but not overly complex in order to avoid over-fitting that leads to high variance error and excessive computing time. Radionuclides releases such as those of ${}^{\mathrm{137}}$Cs from the Fukushima Daiichi NPP are very significant over a long period of time (several weeks) and have a very high temporal variability. In such case, the choice of the discretization is a crucial and challenging task.

In this paper, we explore several ways of improving the source term assessment and its corresponding uncertainty. In the first instance, we address the issue of sampling in the case of massive and highly fluctuating releases of radionuclides to the atmosphere. In the second instance, we propose several statistical models of the likelihood scale matrix depending on the available set of observations.

The outline of the paper is as follows. In Sect. , the observational dataset used in this study and the physical model of the problem are described. We then outline the theoretical aspects of this study and first focus on the concepts of transdimensionality and model selection (Sect. ). The reversible-jump MCMC algorithm (RJ-MCMC) used to shape the release rates vector $\mathit{q}$ is presented in Sect. . Other ways to improve the sampling quality and reduce uncertainties are then explored via the addition of information in Sect. . In particular, methods that combine deposition and concentrations observations or that exploit the temporal and spatial information of the observations in the definition of the likelihood scale matrix $\mathbf{R}$ are proposed. Subsequently, these methods are evaluated in Sect. . Specifically, once the statistical parametrization of the problem is chosen in Sect. , the advantages of the RJ-MCMC are explored in Sect. . The impact of using the deposition measurements and the use of the observation temporal and spatial information are described in Sect.  and . We finally conclude on the contribution of each method.

This article is in line with the authors' earlier studies and is inspired by the insightful work of . In particular, the RJ-MCMC algorithm is applied to the Fukushima Daiichi case by using variable discretizations but a fixed number of steps (i.e. grid cells). Here, it is applied in a more ambitious transdimensional framework, where the number of steps of the discretization is modified through the progression of the RJ-MCMC, and using a significantly larger dataset.

2

The Fukushima Daiichi NPP ${}^{\mathrm{137}}$Cs release and numerical modelling

2.1 Observational dataset

The dataset used in this study consists of ${}^{\mathrm{137}}$Cs air concentration and deposition measurements over the Japanese territory. The vector of ${}^{\mathrm{137}}$Cs air concentration measurements contains $\mathrm{14}\mathrm{248}$ observations from $\mathrm{105}$ stations, whose locations are shown in Fig. . This dataset is significantly larger than the one used by .

Figure 1

Maximum air concentrations of ${}^{\mathrm{137}}$Cs (in Bq m${}^{-\mathrm{3}}$) measured in Japan in March 2011. The light blue points correspond to concentrations below the detection limit. The purple triangle corresponds to the Fukushima Daiichi NPP.

[Figure omitted. See PDF]

This allows us to better evaluate the relevance of our methods and to accurately estimate the uncertainties in the problem. Moreover, the majority of these observations are hourly, which allows for fine sampling of the ${}^{\mathrm{137}}$Cs source term. Additionally, $\mathrm{1507}$ deposition measurements are used. This dataset is described and made available in . The observations are located within $\mathrm{300}\mathrm{km}$ of the Fukushima Daiichi NPP; however, measurements too close to the Fukushima Daiichi plant (typically less than 5 times the spatial resolution of the meteorological fields described in Table ) are not used to consider the short-range dilution effects of the transport model. These measurements represent the averages of each grid cell of a field of deposition observations (in Bq m${}^{-\mathrm{2}}$) measured by aerial radiation monitoring.

2.2 Physical parametrization

The atmospheric ${}^{\mathrm{137}}$Cs plume dispersion is simulated with the Eulerian model ldX, which has already been validated on the Chernobyl and Fukushima Daiichi accidents as well as on the ${}^{\mathrm{106}}$Ru releases in 2017 . The meteorological data used are the 3-hourly high-resolution forecast fields provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) (see Table ).

Table 1

Main set-up parameterizations of the ldX ${}^{\mathrm{137}}$Cs transport simulations. The spatial and temporal resolutions are those of the ECMWF meteorological fields and the transport model. The other entries correspond to parameters of the transport model only.

	ECMWF
Spatial resolution	$\mathrm{0.125}{}^{\circ }\times \mathrm{0.125}{}^{\circ }$
Temporal resolution	3 h
Vertical resolution	11 terrain-following layers (from 0 to 4400 m)
Vertical mixing	K-diffusion parameterized with Louis' closure
	and under unstable planetary boundary layer conditions
Horizontal mixing	Constant horizontal eddy diffusion coefficient, $K_{\mathrm{h}}=\mathrm{0}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
Wet scavenging	$\mathit{\lambda }={\mathrm{\Lambda }}_{\mathrm{0}}{p}_{\mathrm{0}}$ (below-cloud), with ${\mathrm{\Lambda }}_{\mathrm{0}}={\mathrm{5.10}}^{-\mathrm{5}}\mathrm{h}(\mathrm{mm}\mathrm{s}{)}^{-\mathrm{1}}$;
	$\mathit{\lambda }={\mathrm{\Lambda }}_{\mathrm{1}}{p}_{\mathrm{0}}^{\mathrm{0.64}}$ (in-cloud), with ${\mathrm{\Lambda }}_{\mathrm{1}}={\mathrm{5.10}}^{-\mathrm{4}}\mathrm{h}(\mathrm{mm}\mathrm{s}{)}^{-\mathrm{1}}$;
	and $p_{\mathrm{0}}$ is the rainfall intensity
Dry deposition	Constant deposition velocity, $v_{\mathrm{d}}={\mathrm{2.10}}^{-\mathrm{3}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$

The radionuclides, in particular ${}^{\mathrm{137}}$Cs, were mainly deposited in the Fukushima Prefecture, north-west of the Fukushima Daiichi NPP. The bulk of the releases occurred over a 2-week period beginning on 11 March. Therefore, the model simulations start on 11 March 2011 at 00:00 UTC and end on 24 March 2011, corresponding to 312 h during which radionuclides could have been released. The release is assumed to be spread over the first two vertical layers of the model (between 0 and 160 m). The height of the release has a small impact, as the modelled predictions are compared with observations sufficiently far away from the Fukushima Daiichi NPP and the release mostly remains within the atmospheric boundary layer. Finally, it is assumed that all parameters describing the source aside from the reconstructed ones are known. This is, for instance, the case for the release height.

3 Methodological aspects of the inverse problem

3.1 Reversible-jump MCMC

In this section, we describe the Bayesian reversible-jump MCMC algorithm used in the following to reconstruct the vector $\mathit{q}$ discretized over a variable adaptive grid. A good description of the RJ-MCMC, which was introduced by , can be found in . The algorithm is also clearly explained and applied by . In particular, the method was used by in the field of inverse modelling for the assessment of substance releases as well as by to quantify Fukushima Daiichi releases, although with a fixed number of variable-sized grid cells. A key asset of the RJ-MCMC technique is its ability to sample a highly variable release by balancing bias and variance errors.

The RJ-MCMC algorithm is a natural extension of the traditional Metropolis–Hastings (MH) algorithm to transdimensional discretization grids. In a traditional MH, the distribution of the release is assessed via a succession of random walks e.g.. In particular, special random walks allow the pdf of the discretization steps within the control vector $\mathit{x}$ to be sampled.

Indeed, as described above, the function that describes the evolution of the release rate over time is a vector or step function, where each step (i.e. a grid cell) defines a constant release rate over a time interval. When using the RJ-MCMC, these time intervals are neither fixed nor regularly spaced. Hence, an objective is to retrieve the best stepwise time partition of the source term, which defines the release rates.

To do this, we use the concept of boundary: a boundary is a time that separates two constant release rates (i.e. that separates two grid cells where the source term is defined).

Figure 2

Example of a source term partition: $\mathbf{\Lambda }=({\mathit{\lambda }}_{\mathrm{0}},\mathrm{\dots },{\mathit{\lambda }}_{\mathrm{6}}$) is the set of boundaries, and $\mathit{q}=(q_{\mathrm{0}},\mathrm{\dots },q_{\mathrm{5}})$ is the vector of release rates. The release rate defined on the time interval $[{\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{3}}]$ is constant and equal to $q_{\mathrm{2}}$.

[Figure omitted. See PDF]

Let $\mathbf{\Lambda }$ be the set of $N_{\mathrm{b}}$ boundaries. Boundaries are $\mathrm{0}$, $\mathrm{312}$, and integers in $\mathit{\{}\mathrm{1},\mathrm{2}\mathrm{\dots },\mathrm{311}\mathit{\}}$ ($\mathrm{0}$ and $\mathrm{312}$ are natural and fixed boundaries), i.e. the boundaries are a selection among the hours of the release. Finding the probability distribution of the time steps associated with $\mathit{q}$ is equivalent to finding the probability distribution of $\mathbf{\Lambda }$. To retrieve the best discretizations of $\mathit{q}$, we can therefore sample $\mathbf{\Lambda }$. Figure  presents an example of a source term partition.

The RJ-MCMC is basically equivalent to an MH algorithm except that some iterations of the algorithm are not intra-model random jumps (i.e. at constant $\mathbf{\Lambda }$) but rather possibly inter-model random jumps (i.e. at non-constant $\mathbf{\Lambda }$) or transdimensional (i.e. non-constant $\mathbf{\Lambda }$ size).

Specifically, we use three types of random jumps that govern the positions and number of boundaries: birth, death, and move. For these three processes, we draw on the work of , , , and (importantly) . The birth process is the creation of a new boundary not yet present in $\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{311}\mathit{\}}$, the death process is the removal of an existing boundary in $\mathbf{\Lambda }$, and the move process is the displacement of an existing boundary in $\mathbf{\Lambda }$. These processes must respect the detailed balance criterion in order to ensure the convergence of the Markov chain to the target distribution .

3.2 Handling observations and their errors

3.2.1 Measurements fusion

In this subsection, we propose a method to factor both concentration and deposition measurements into a Bayesian sampling by simply taking their disparity in the observation error matrix $\mathbf{R}$ into account. First of all, $\mathbf{R}$ is diagonal and is, therefore, described by a collection of scale parameters. Specifically, two scale parameters, $r_{\mathrm{c}}$ and $r_{\mathrm{d}}$, are used when modelling $\mathbf{R}$ and are associated with the air concentration and deposition measurements respectively. Thus we can write the following:

2$$\mathbf{R}=\mathrm{diag}(r_{\mathrm{c}},\mathrm{\dots },r_{\mathrm{c}},r_{\mathrm{d}},\mathrm{\dots },r_{\mathrm{d}}),$$ where $\mathrm{diag}\left(\mathit{\right\{}{r}_i{\mathit{\}}}_i)$ is the diagonal matrix of diagonal entries $\mathit{\{}{r}_i{\mathit{\}}}_i$. Note that we use the sorting algorithm described by , which is applied to both the air concentration and deposition measurements, in all subsequent applications. This algorithm assigns a different scale parameter $r_{\mathrm{np}}$ to certain non-pertinent observations to prevent them from artificially reducing $r_{\mathrm{c}}$ and $r_{\mathrm{d}}$. With this definition of $\mathbf{R}$, the scale variables $r_{\mathrm{c}},r_{\mathrm{d}}$, and $r_{\mathrm{np}}$ are included in the source vector $\mathit{x}$ and, therefore, sampled using a MCMC.

3.2.2 Spatio-temporal observation error scale matrix modelling

In this subsection, we propose modelling $\mathbf{R}$ using the spatio-temporal information of the concentration observations. We assume that the error scale parameter corresponding to an observation–prediction pair depends mainly on the error in the model predictions. We also propose considering a prediction (i.e. an estimation of the number of radionuclides present at a definite time and place) as a point in a radionuclide plume at a specific time. The prediction modelling error varies according to its spatio-temporal position. Spatially, if two points of interest are distant by ${\mathrm{\Delta }}_x$, the difference in the errors in their attached model prediction is empirically assumed to be proportional to ${\mathrm{\Delta }}_x$. Temporally, let us consider two points in a plume with coordinates $x_{\mathrm{1}},x_{\mathrm{2}}$ at times $t$ and $t+{\mathrm{\Delta }}_t$ respectively, where $x_{\mathrm{1}}$ and $x_{\mathrm{2}}$ are the longitude and latitude. These two points are moving with the plume, and they are spatially distant by $v\times {\mathrm{\Delta }}_t$ at time $t$, where $v$ is a reference wind speed, representing the average wind over the accidental period. Therefore, the difference in the modelling errors of these points is estimated to be proportional to ${\mathrm{\Delta }}_x=v\times {\mathrm{\Delta }}_t$. In that respect, for each concentration observation, $y_i$ can be associated with three coordinates:

3$$\left(x_{\mathrm{1},i},x_{\mathrm{2},i},v\cdot \left(t_i+{\displaystyle \frac{{\mathrm{\Delta }}_{t,i}}{\mathrm{2}}}\right)\right)$$ where $x_{\mathrm{1},i}$, $x_{\mathrm{2},i}$ are the coordinates (in km) of the observation, $v$ is set equal to $\mathrm{12}\mathrm{km}{\mathrm{h}}^{-\mathrm{1}}$, and $t_i$ and ${\mathrm{\Delta }}_{t,i}$ are the starting time and the duration of the observation in hours respectively. This average wind speed value has not been extensively researched and is only used to estimate the potential of the method. A clustering algorithm such as the $k$-means algorithm can be applied to partition these three coordinates' observation characteristics into several groups. A scale parameter $r_{\mathrm{c},i}$ can then be assigned to each group. If the algorithm for sorting pertinent and non-pertinent observations is applied, this yields 4$$\mathbf{R}=\mathrm{diag}(r_{\mathrm{c},i},r_{\mathrm{np}},r_{\mathrm{c},j},\mathrm{\dots },r_{\mathrm{d}},r_{\mathrm{np}}),$$ where $r_{\mathrm{c},i},r_{\mathrm{d}}$, and $r_{\mathrm{np}}$ are the scale parameters associated with the spatial or spatio-temporal clusters of air concentration observations, the deposition measurements, and the non-pertinent observations respectively. The case where only the spatial qualification of the observation is used instead of both the spatial and temporal qualification is also addressed in this study.

4

Application to the Fukushima Daiichi NPP ${}^{\mathrm{137}}$Cs release

In this section, we apply the previously introduced methods to the reconstruction of the Fukushima Daiichi ${}^{\mathrm{137}}$Cs source term between 11 and 24 March. However, we first define the key statistical assumptions and parameters of the inversion set-up.

4.1 Statistical parametrization

The prior distributions determining the posterior distribution and the transition probabilities of the RJ-MCMC algorithm are defined in the next two sections.

4.1.1 Definition of the prior density functions

In this section, we specify the prior distributions of the source vector components: release rates and hyperparameters. The source release rate logarithms are assumed to follow a folded Gaussian prior distribution, so as to concurrently constrain the number of boundaries and the magnitude of the release rates. In the general case with $N_{\mathrm{imp}}=\mathrm{312}$, the scale matrix $\mathbf{B}$ of the folded Gaussian prior is defined as a $\mathrm{312}\times \mathrm{312}$ diagonal matrix. The entry $b_k$ of $\mathbf{B}$ is associated with the $k$th hour in $\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{311}\mathit{\}}$, comprising hours where we consider the release of ${}^{\mathrm{137}}$Cs possible. This yields the following:

5$$p\left(\ln \mathit{q}\right|N_{\mathrm{imp}})=\prod _{k=\mathrm{0}}^{N_{\mathrm{imp}-\mathrm{1}}}\sqrt{{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }{b}_k}}}\left(e^{-\frac{(\ln q_k-\ln q_{\mathrm{b}}{)}^{\mathrm{2}}}{\mathrm{2}{b}_k}}\right),$$ where the folding position is defined as equal to the location term $\ln q_{\mathrm{b}}$. We divide the 312 time intervals into two groups. The first “unconstrained” group consists of the time intervals during which we have noticed that a radionuclide emission does not trigger any noticeable observation. In other words, it is a time interval during which we have no information on a potential release. The unconstrained time intervals are the first 17 h of 11 March (hours in $\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{16}\mathit{\}}$), the days of 16 and 17 March (hours in $\mathit{\{}\mathrm{116},\mathrm{117},\mathrm{\dots },\mathrm{165}\mathit{\}}$), and the last 18 h of 23 March (hours in $\mathit{\{}\mathrm{294},\mathrm{295},\mathrm{\dots },\mathrm{311}\mathit{\}}$). The second group consists of the constrained time intervals. Two hyperparameters – $b_{\mathrm{c}}$ for the constrained time intervals and $b_{\mathrm{nc}}$ for the unconstrained time intervals – are used to describe the matrix $\mathbf{B}$. If, instead of two, a single hyperparameter $b$ is used for all release rates, the unconstrained time intervals are neither regularized by the observations nor by the prior. Indeed, to account for high release rates, samples of the distribution of $b$ reach significant values. Low-magnitude release rates are then constrained very little by the corresponding prior. Release rates that are neither constrained by the observations nor by the priors can, in turn, reach very large values. This phenomenon was put forward by . The mean of the folded Gaussian prior is chosen as $\ln q_{\mathrm{b}}=\mathrm{0}$ in order to regularize the unconstrained intervals. It is important to note that a release rate $\ln q_i$ sampled by the RJ-MCMC is an aggregate of these hourly release rates, with each of them associated with either $b_{\mathrm{c}}$ or $b_{\mathrm{nc}}$. Hence, we define the scale of the prior associated with the release rate $\ln q_i$ defined between times $t_i$ and $t_i+k_i{\mathrm{\Delta }}_t$ as follows: 6$$\begin{array}{rl}{b}_i& ={\displaystyle \frac{{\sum }_{j=t_i}^{t_i+k_i{\mathrm{\Delta }}_t}{\mathbf{B}}_{j,j}}{k_i{\mathrm{\Delta }}_t}}={\displaystyle \frac{n_{\mathrm{c},i}{b}_{\mathrm{c}}+(k_i{\mathrm{\Delta }}_t-n_{\mathrm{c},i})b_{\mathrm{nc}}}{k_i{\mathrm{\Delta }}_t}}\\ & =w_{\mathrm{c},i}{b}_{\mathrm{c}}+w_{\mathrm{nc},i}{b}_{\mathrm{nc}},\end{array}$$ where ${\mathbf{B}}_{j,j}$ is the $j$th diagonal coefficient of $\mathbf{B}$, $n_{\mathrm{c},i}$ is the portion of constrained time intervals in $k_i{\mathrm{\Delta }}_t$, and $w_{\mathrm{c}{,}_i}$ and $w_{\mathrm{nc},i}$ are the weighted numbers of the respective constrained and unconstrained hourly release rates included in the release rate $\ln q_i$. This definition is reflected in the prior cost: $$\begin{array}{}\text{7}& {\displaystyle }{\mathbf{J}}_{\mathrm{prior},\ln \mathit{q}}=\sum _{i=\mathrm{0}}^{N_{\mathrm{imp}}-\mathrm{1}}\left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln \left({\displaystyle \frac{\mathit{\pi }{b}_i}{\mathrm{2}}}\right)+{\displaystyle \frac{{\left(\ln q_i\right)}^{\mathrm{2}}}{\mathrm{2}{b}_i}}\right\}\text{8}& {\displaystyle }\begin{array}{rl}=& \sum _{i=\mathrm{0}}^{N_{\mathrm{imp}}-\mathrm{1}}\left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln \left({\displaystyle \frac{\mathit{\pi }(w_{\mathrm{c},i}{b}_{\mathrm{c}}+w_{\mathrm{nc},i}{b}_{\mathrm{nc}})}{\mathrm{2}}}\right)\right.\\ & \left.+{\displaystyle \frac{{\left(\ln q_i\right)}^{\mathrm{2}}}{\mathrm{2}(w_{\mathrm{c},i}{b}_{\mathrm{c}}+w_{\mathrm{nc},i}{b}_{\mathrm{nc}})}}\right\},\end{array}\end{array}$$ where $N_{\mathrm{imp}}=N_{\mathrm{b}}-\mathrm{1}$ is the number of time steps ($N_{\mathrm{b}}$ being the number of boundaries) considered at a certain iteration by the RJ-MCMC and which characterizes the grid on which $\ln \mathit{q}$ is defined. The prior scale parameter $b_{\mathrm{c}}$ is included in the source vector $\mathit{x}$ and sampled with the rest of the variables. Such a prior on the release rate logarithms also has a constraining effect on the model's complexity (Occam's razor principle) through the normalization constant of the prior pdf.

An exponential prior distribution is used for the boundaries $\mathbf{\Lambda }=({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_k)$: 9$$\begin{array}{rl}p& ({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_k)=\\ & \left\{\begin{array}{ll}\frac{e^{-k}}{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{b},\max }}\frac{N_{\mathrm{b},\max }\mathrm{!}}{i\mathrm{!}(N_{\mathrm{b},\max }-i)\mathrm{!}}{e}^{-i}}& \text{if}k\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots };N_{\mathrm{b},\max }\mathit{\}};\\ \mathrm{0}& \text{otherwise}.\end{array}\right.\end{array}$$ Here, $N_{\mathrm{b},\max }=\mathrm{311}$. This prior has the effect of penalizing models that are too complex.

Furthermore, because no a priori information is available, the prior distributions on the coefficients of the observation error scale matrix are assumed to be uniform, and the prior values on the regularization scale terms are assumed to follow Jeffreys' prior distribution . The lower and upper bounds of the uniform prior on the coefficients of the observation error matrix are set as zero and a large value, the latter of which is not necessarily realistic, as in .

4.1.2 Parameters of the MCMC algorithm

In the Markov chain, the variables describing the source $\mathit{x}$ are initialized randomly. The intra-model (i.e. with a fixed boundary partition) transitions defining the stochastic walks are set independently for each variable. Each transition probability is defined as a folded normal distribution, following . Locations of the transition probabilities are the values of the variables at the current step. Variances of the related normal distribution are initialized at ${\mathit{\sigma }}_{\ln q}=\mathrm{0.3}$ and ${\mathit{\sigma }}_r=\mathrm{0.01}$, with the values being chosen by means of experiments. However, the value of ${\mathit{\sigma }}_r$ is adaptive and is updated every $\mathrm{1000}$ iterations according to the predictions' values. Furthermore, the value of ${\mathit{\sigma }}_r$ is not uniform across groups of observations: there are as many variances as observation clusters.

The algorithm runs initially as a classic MH for ${\mathrm{10}}^{\mathrm{4}}$ iterations (only intra-model walks are allowed at this stage). Indeed, allowing transdimensional walks at the beginning of the run sometimes cause Markov chains to get stuck in local minima. In total, $\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ iterations are used; this large number ensures algorithm convergence (i.e. sufficient sampling of the posterior distribution of $\mathit{x}$) and corresponds to approximately 1 d of calculation for a $\mathrm{12}$-core computer. The burn-in is set to ${\mathrm{10}}^{\mathrm{6}}$ iterations.

The inter-model transdimensional walks (with changing boundaries) are described in Appendix , as they are technical. For this transdimensional case, $u_{\ln q}$ is the same for birth and death processes and is equal to $\mathrm{0.3}$. There is no variance parameter for the move walk. The probability of proposing an intra-model transition, a birth transition, a death transition and a move transition is $\mathrm{1}/\mathrm{2}$, $\mathrm{1}/\mathrm{6}$, $\mathrm{1}/\mathrm{6}$, and $\mathrm{1}/\mathrm{6}$ respectively.

4.2 Results

To evaluate the quality of our reconstructed source term as well as the impact of the techniques proposed in Sect. , the distributions of the ${}^{\mathrm{137}}$Cs source are sampled and compared using several settings:

• Sect.  provides an activity-concentration-based reconstruction of the ${}^{\mathrm{137}}$Cs source with a log-Cauchy likelihood and a threshold (to ensure it is defined for zero observations or predictions) of $\mathrm{0.5}$ Bq m${}^{-\mathrm{3}}$ , the application of the observation sorting algorithm, and a spatial representation of $\mathbf{R}$ with 10 parameters. The vector

  10$$\mathit{x}=\left(\ln \mathit{q},r_{\mathrm{c},\mathrm{1}},\mathrm{\dots },r_{\mathrm{c},\mathrm{9}},r_{\mathrm{np}},b_{\mathrm{c}},({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N_{\mathrm{imp}}-\mathrm{1}})\right)$$is sampled. A log-Cauchy likelihood is chosen because the associated reconstruction produces the best model-to-data comparisons among diverse choices of likelihoods and representations of $\mathbf{R}$. The sampled source is compared to the source terms of and , and all variable marginal distributions are described.

• In Sect. , we are interested in quantifying the impact of the addition of deposition measurements in the sampling with the method presented in Sect. . The source term is sampled in two configurations, namely with or without the addition of deposition measurements. That is, the vector11$$\begin{array}{rl}\mathit{x}& =\left(\ln \mathit{q},r_{\mathrm{c},\mathrm{1}},\mathrm{\dots },r_{\mathrm{c},\mathrm{9}},r_{\mathrm{np}},r_{\mathrm{d}},b_{\mathrm{c}},({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },\right.\\ & \left.{\mathit{\lambda }}_{N_{\mathrm{imp}}-\mathrm{1}})\right)\end{array}$$is sampled with or without $r_{\mathrm{d}}$. A threshold of $\mathrm{0.5}$ Bq m${}^{-\mathrm{2}}$ is chosen for the deposition measurements.

• Finally, we draw a selection of violin plots in Sect.  to assess the impact of the distribution of the observations in space and time clusters using the approach introduced in Sect. . More precisely, the fit to the observations of the reconstructed source is evaluated given the type of clustering (spatial or spatio-temporal) and the number of clusters.

4.2.1 Best reconstruction and comparison to other source terms

In this section, we provide the most appropriate reconstruction of the source term, based on the activity concentration measurements, from several choices of likelihoods and $\mathbf{R}$ representations. This best reconstruction is chosen with the help of a model-to-data comparison, i.e. FAC (factor) scores, representing the proportion of measurements for which the observed and modelled values agree. The chosen likelihood is a log-Cauchy distribution with a threshold of $\mathrm{0.5}$ Bq m${}^{-\mathrm{3}}$, and air concentration measurements are spatially clustered in nine groups using the $k$-means method.

Figure  shows the temporal evolution of the median release rate with its associated standard deviation (i.e. the evolution of the hourly release rates samples that are obtained from the $N_{\mathrm{imp}}$ RJ-MCMC samples). This temporal evolution is compared to the source term from , which was estimated from gamma dose rate measurements but from the same meteorological data and transport model, and to the more recent source term from .

Figure 3

Evolution of the release rate (in Bq s${}^{-\mathrm{1}}$) describing the ${}^{\mathrm{137}}$Cs source (UTC) sampled with the RJ-MCMC algorithm. The blue line corresponds to the sampled median release rate. The light blue area corresponds to the area between ${\mathit{\mu }}_{\ln q}$, the median, and ${\mathit{\sigma }}_{\ln q}$, the standard deviation, of the reconstructed hourly release rates: $\left[{\mathit{\mu }}_{\ln q}-{\mathit{\sigma }}_{\ln q},{\mathit{\mu }}_{\ln q}+{\mathit{\sigma }}_{\ln q}\right]$. Our source term is compared to the source terms of and .

[Figure omitted. See PDF]

Firstly, we observe that a good match is obtained between the median release and the source term of or (on the constrained time intervals), as the differences between releases under ${\mathrm{10}}^{\mathrm{9}}$ Bq s${}^{-\mathrm{1}}$ can be neglected. Indeed, between the 11 and the 19 March 2011, the three source terms have similar peak magnitudes (close to ${\mathrm{10}}^{\mathrm{11}}$ Bq s${}^{-\mathrm{1}}$); these peaks are observed at neighbouring time intervals. Between the 19 and the 21 March, a notable difference can be observed: the RJ-MCMC source term predicts very high 1 h peaks (close to ${\mathrm{10}}^{\mathrm{12}}$ Bq s${}^{-\mathrm{1}}$), whereas the source term distributes the release uniformly over the whole period. Secondly, release rates are subject to large variations. Indeed, there are periods of important temporal variability in the release rate (e.g. between 19 and 21 March) and of low temporal variability (e.g. between 11 and 14 March). This justifies the use of the RJ-MCMC transdimensional algorithm, as it allows one to reconstruct certain time intervals of the release with a greater accuracy and others with less complexity. Thirdly, the temporal evolution of the release can be explored, and several prominent peaks can be observed with a magnitude approaching or surpassing ${\mathrm{10}}^{\mathrm{11}}$ Bq s${}^{-\mathrm{1}}$: between 12 and 14 March, between 14 and 16 March, on 18 March, and between 19 and 21 March.

Figure 4

Densities or averages of variables describing the ${}^{\mathrm{137}}$Cs source reconstructed with the RJ-MCMC technique: (a) density of relative differences in log${}_{\mathrm{10}}$ between observations and predictions; (b) density of the total retrieved released activity (TRRA) of ${}^{\mathrm{137}}$Cs during the Fukushima Daiichi accident (in Bq) – outliers are removed; (c) histogram of the number of boundaries; and (d) mean ($\pm$ standard deviation) of the number of boundaries as a function of time (UTC).

[Figure omitted. See PDF]

Figure a describes the histogram of relative differences in log${}_{\mathrm{10}}$ between predictions (medians of the predictions computed using the samples are employed) and air concentration observations larger than $\mathrm{0.5}$ Bq m${}^{-\mathrm{3}}$. A value of $-\mathrm{1}$, therefore, corresponds to an observation 10 times larger than a prediction. A good fit between observations and predictions can be observed, although the observations are globally underestimated by the predictions.

Figure b shows the pdf of the total release (in Bq). It can be seen from this density that most of the total release distribution mass is sampled between $\mathrm{10}$ and $\mathrm{20}$ PBq and peaks at $\mathrm{14}$ PBq. This is consistent with previous studies which estimate the release to be between $\mathrm{5}$ and $\mathrm{30}$ PBq . However, note that these works attempted to estimate potential releases beyond 24 March.

Figure c plots the pdf of the number of boundaries (i.e. a measure of model complexity) sampled by the algorithm. Here, the minimum of the distribution is 25 boundaries and the maximum is 45 with a peak at 35 boundaries. Recall that the maximum number of boundaries is 313.

Figure d shows the average number of sampled boundaries at each time and the associated standard deviation. The $y$ axis corresponds to the number of boundaries counted: for a given hour $t$, we draw 12$$E\left[\mathrm{Card}\mathit{\right\{}{\mathit{\lambda }}_i\in [t-\mathrm{5}\mathrm{\dots }t+\mathrm{5}]|{\mathit{\lambda }}_i\in ({\mathit{\lambda }}_{\mathrm{0}},\mathrm{\dots },{\mathit{\lambda }}_{N_{\mathrm{imp}}}\left)\mathit{\right\}}],$$ i.e. the average over all samples of the number of boundaries counted around $t$. In essence, the curve shows the number of boundaries necessary to model the release at a certain time (i.e. the variability in the release according to the algorithm). We observe that there is a large variability between 14 and 15 March as well as between 19 and 21 March. This correlates very well with the observable results in Fig. . This observation once again confirms the relevance of using the RJ-MCMC algorithm: periods of high variability can be sampled finely.

4.2.2 Impact of using deposition measurements

In this section, we study the impact of adding $\mathrm{1507}$ deposition measurements to the $\mathrm{14}\mathrm{248}$ air concentration measurements. A log-Cauchy distribution with a threshold of $\mathrm{0.5}$ Bq m${}^{-\mathrm{3}}$ for the air concentration measurements and $\mathrm{0.5}$ Bq m${}^{-\mathrm{2}}$ for the deposit measurements is employed. Air concentration measurements are clustered with a spatial $k$-means algorithm in nine groups, as described in Sect. . To quantify the impact of the added measurements, we sample the following vector with an RJ-MCMC algorithm in two cases, with or without deposition measurements, i.e. with or without $r_{\mathrm{d}}$, as follows:

13$$\mathit{x}=\left(\ln \mathit{q},r_{\mathrm{c},\mathrm{1}},\mathrm{\dots },r_{\mathrm{c},\mathrm{9}},r_{\mathrm{np}},r_{\mathrm{d}},b_{\mathrm{c}},({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N_{\mathrm{imp}}-\mathrm{1}})\right).$$ In the case without deposition measurements, we use the same samples as in Sect. . A posteriori marginals of these control variables are displayed in Fig. .

Figure 5

The panels correspond to probability densities or averages of variables describing the ${}^{\mathrm{137}}$Cs source reconstructed with the RJ-MCMC technique using air concentration and deposition measurements. Panel (a) displays the medians and standard deviations of hourly release rates. The red and blue curves correspond to the sampled median release rates. The red and light blue areas correspond to the areas between ${\mathit{\mu }}_{\ln q}$, the median, and ${\mathit{\sigma }}_{\ln q}$, the standard deviation, of the sampled hourly release rates for each reconstruction: $\left[{\mathit{\mu }}_{\ln q}-{\mathit{\sigma }}_{\ln q},{\mathit{\mu }}_{\ln q}+{\mathit{\sigma }}_{\ln q}\right]$. Panel (b) shows the probability density of the total release of ${}^{\mathrm{137}}$Cs during the Fukushima Daiichi accident (in Bq). Panel (c) provides a histogram of the number of boundaries. Panel (d) shows the mean ($\pm$ standard deviation) of the number of boundaries according to the hour considered.

[Figure omitted. See PDF]

In Fig. a, the medians of the two sampling datasets are similar, except for the release rates between 11 and 14 March. Indeed, when using air concentration and deposition measurements, the release rate peaks on 12 March instead of on 13 March, the latter of which is the case when using air concentrations alone.

In Fig. b, the main mode of the total retrieved released activity (TRRA) distribution retrieved using the air concentration and deposition measurements is $\mathrm{5}$ PBq larger than the TRRA reconstructed with the air concentration measurements alone. Moreover, the reconstructed density with both types of measurements has a large variance and leads to the possibility of a larger release (up to $\mathrm{75}$ PBq). The density of the number of boundaries with both types of measurements is shifted to higher values compared with the density with air concentrations alone, which indicates a reconstructed release of greater complexity (Fig. c). This is consistent with the fact that the second release is reconstructed with a larger number of observations.

In Fig. d, the number of boundaries was estimated to be half as large on 14 March when deposition measurements were used, but it was estimated to be larger on 15 March. These differences are rather surprising insofar as there is no obvious disparity between the medians of the release rates for the two reconstructions in Fig. a on these days.

In Table , we evaluate the FAC$$2, FAC$$5, and FAC$$10 scores for both cases. FAC scores are common metrics for measuring the fit of a source term with available observations. The FAC$$2, FAC$$5, and FAC$$10 are the proportion of measurements for which the observed and modelled values agree by a factor of 2, 5, or 10 respectively. For example, a FAC$$2 score of $\mathrm{0.32}$ with respect to the air concentration measurements $(y_i{)}_{i\in \mathit{\{}\mathrm{1},N_{\mathrm{obs}}\mathit{\}}}$ means that $\mathrm{32}\mathit{\%}$ of the predictions $\left(\right(\mathbf{H}\mathit{q}{)}_i{)}_{i\in \mathit{\{}\mathrm{1},N_{\mathrm{obs}}\mathit{\}}}$ (i.e. simulated measurements) are such that 14$$\mathrm{0.5}<{\displaystyle \frac{(\mathbf{H}\mathit{q}{)}_i+y_{\mathrm{t}}}{y_i+y_{\mathrm{t}}}}<\mathrm{2},$$ where $N_{\mathrm{obs}}$ is the number of observations considered, and $y_{\mathrm{t}}$ is a threshold set at $\mathrm{0.5}$ Bq m${}^{-\mathrm{3}}$. This threshold is added to limit the weight of low values, considered of lesser interest, in the indicator. Finally, to calculate the FAC scores, only observations above $\mathrm{1}$ Bq m${}^{-\mathrm{3}}$ are considered.

Table 2

Comparison of observed and simulated measurements from reconstructed predictions with air concentration measurements only or with air concentration and deposition measurements. The FAC$$2, FAC$$5, and FAC$$10 are the proportion of measurements for which the observed and modelled values agree by a factor of 2, 5, or 10 respectively. For example, the table reads as follows: the predictions calculated based on air concentrations only obtain a FAC$$5 score of 0.58 with respect to the deposition measurements.

Origin of the source term	FAC$$2	FAC$$5	FAC$$10
RJ-MCMC (based on air concentrations only)
Air concentration	0.32	0.68	0.85
Deposition	0.31	0.58	0.75
RJ-MCMC (air concentration and deposition)
Air concentration	0.31	0.67	0.85
Deposition	0.49	0.77	0.86
Source term of
Air concentration	0.18	0.52	0.76
Deposition	0.28	0.64	0.74
Source term of
Air concentration	0.17	0.52	0.78
Deposition	0.29	0.64	0.76

The FAC scores show that the source term obtained with the RJ-MCMC has a better agreement with the air concentration measurements than the or source terms. This shows that the RJ-MCMC is able to sample source terms consistent with air concentration measurements but does not allow us to reach a conclusion regarding its performance: the sampling algorithm was optimized using air concentrations, whereas, for example, the source term was recovered from dose rate measurements. However, the source term found with the RJ-MCMC on air concentration measurements can be fairly evaluated using deposition measurements. It is observed that the fits with the deposition measurements of the three source terms are similar.

Furthermore, as expected, Table 2 shows that there is a significant difference in the fit with the deposition measurements when they are added to the observation dataset used to reconstruct the source term. However, the addition of the deposition measurements has no impact on the quality of the fit to the air concentration measurements. This suggests that an RJ-MCMC optimized on dose rate measurements, in addition to the deposition and air concentrations, may obtain a good agreement with the dose rate measurements, without degrading its results on deposition or air concentrations.

4.2.3

Impact of the observation error matrix $\mathbf{R}$ representation

In this section, we investigate the spatial and spatio-temporal representation of the $\mathbf{R}$ likelihood scale matrix. Figure  presents violin plots of the FAC$$2 and FAC$$5 scores for air concentration measurements against the number of scale parameters chosen to describe $\mathbf{R}$ (i.e. the number of observations clusters). Violin plots are created out of the best scores of 44 RJ-MCMC sampling datasets with diverse likelihoods chosen after in order to make the conclusions less dependent on statistical modelling. The chosen numbers of spatial clusters (1, 2, and 9) and spatio-temporal clusters (3 and 7) of the violins are a selection for the sake of clarity.

Figure 6

Violin plots of the air concentration FAC$$2 and FAC$$5 scores against the observation error scale parameters' spatio-temporal groups number. Points composing the violin plots consist of the best FAC$$2 and FAC$$5 scores of 44 sampling datasets using diverse configurations of likelihoods and thresholds (log-normal, log-Laplace, and log-Cauchy) chosen according to . The horizontal segments within the violins represent (from bottom to top) the respective first, second, and third quartiles. The chosen numbers of spatial and spatio-temporal clusters are only a selection for the sake of clarity.

[Figure omitted. See PDF]

For both the FAC$$2 and FAC$$5 air concentration scores, it can be observed that the spatio-temporal labelling of the observations does not increase the quality of the scores, whereas the number of spatial clusters does. For example, we can observe a difference of $\mathrm{0.02}$ between the second quartile of the FAC$$2 score for sampling datasets with two spatial groups and the second quartile of the FAC$$2 score for sampling datasets with nine spatial groups. One hypothesis is that the error is mainly related to the distance to the source and remains relatively homogeneous in time. However, the use of an accurate estimate of the effective wind speed as a function of time and depending on the position of the observations is necessary to draw conclusions; such a study is outside the scope of this paper.

5 Summary and conclusions

In this paper, we investigated a transdimensional sampling method to reconstruct highly fluctuating radionuclide atmospheric sources and applied it to assess the ${}^{\mathrm{137}}$Cs Fukushima Daiichi release. Furthermore, we proposed two more methods to add information to the model in order to reduce uncertainties attached to such a complex release.

Firstly, an RJ-MCMC algorithm was defined. This MH extension to transdimensional meshes allows one to sample both the source term and its underlying adaptive discretization in time. In the case of complex releases, an adaptive grid allows one to reconstruct the source term with high accuracy and discretization consistent with the available observations.

Secondly, we have focused on ways to add information by adapting the representation of the observation error matrix $\mathbf{R}$. We proposed the assimilation of deposition measurements or spatial and temporal information on the air concentration observations by increasing the complexity of $\mathbf{R}$.

Subsequently, the distributions of the variables defining the source of ${}^{\mathrm{137}}$Cs were fully sampled. This enabled the estimation of the uncertainties associated with these variables as well as the evaluation and demonstration of the merits of the methods.

Firstly, we observed that the predictions reproduce the air concentration and deposition observations well. When using air concentrations only, the FAC$$2 and FAC$$5 scores corresponding to air concentration measurements are $\mathrm{0.32}$ and $\mathrm{0.68}$ respectively. The main part of the best total release distribution mass is estimated to be between $\mathrm{10}$ and $\mathrm{20}$ PBq and peaks at $\mathrm{14}$ PBq, which is in accord with previous estimations.

Secondly, the periods of high variability in ${}^{\mathrm{137}}$Cs releases have been reconstructed with accuracy by the RJ-MCMC algorithm. The transdimensional method also allowed the periods of low temporal variability to be sampled more coarsely, thereby avoiding variance errors and saving computing time, which, in particular, confirms the conclusions of . The priors that we chose allowed for an efficient regularization of the release rates and their adaptive mesh.

Finally, the use of combined air concentration and deposition measurements had an impact on the reconstruction of the TRRA distribution, the mean of which increased by $\mathrm{10}$ PBq. The spatial clustering methods proved to increase the quality of the predictions by improving the model-to-data comparison.

We recommend the use of the RJ-MCMC for long-release source reconstruction. It allows one to sample highly variable source terms with great accuracy by solving the bias–variance trade-off. We also recommend the use of advanced modelling of $\mathbf{R}$ considering the deposition measurements and the spatial information of the air concentration observations which have proven to be valuable in reducing the observation–prediction discrepancies.

Appendix A RJ-MCMC mathematical details

Here, we describe the mathematical details of the three inter-model walks: move, birth, and death. The move process is intra-dimensional. Provided that the move is a symmetric motion, the associated $g$-transition distribution is symmetric, i.e. A1$${\displaystyle \frac{g\left({\mathit{x}}_i\right|{\mathit{x}}_j)}{g\left({\mathit{x}}_j\right|{\mathit{x}}_i)}}=\mathrm{1}.$$

The birth process is the creation of a new boundary, a walk from ${\mathit{x}}_i$ (a vector composed of $n$ boundaries) to ${\mathit{x}}_j$ (a vector composed of $n+\mathrm{1}$ boundaries) is proposed. In other words, by generating a new boundary, one release rate is destroyed and two new release rates emerge from this destruction. We therefore assign two new release rates $\ln q_{k-\mathrm{1}}^{\prime }=\ln q_{k-\mathrm{1}}-u_{\ln q}$ and $\ln q_k^{\prime }=\ln q_{k-\mathrm{1}}+u_{\ln q}$, with $u_{\ln q}$ (a Gaussian noise). The transition probability is defined as follows : A2$$\begin{array}{rl}g& ({\mathbf{\Lambda }}_i,\ln {\mathit{q}}_i|{\mathbf{\Lambda }}_j,\ln {\mathit{q}}_j)=g({\mathbf{\Lambda }}_i|{\mathbf{\Lambda }}_j,\ln {\mathit{q}}_j)\\ & g\left(\ln {\mathit{q}}_i\right|{\mathbf{\Lambda }}_i,{\mathbf{\Lambda }}_j,\ln {\mathit{q}}_j)=g({\mathbf{\Lambda }}_i\left|{\mathbf{\Lambda }}_j\right)\\ & g\left(\ln {\mathit{q}}_i\right|\ln {\mathit{q}}_j,{\mathbf{\Lambda }}_i).\end{array}$$ The probability of a birth at a certain position among $N_{\mathrm{b},\max }-N_{\mathrm{imp}}+\mathrm{1}$ positions not occupied by boundaries is A3$$g\left({\mathbf{\Lambda }}_i\right|{\mathbf{\Lambda }}_j)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{b},\max }-N_{\mathrm{imp}}+\mathrm{1}}}.$$ Furthermore, the probability of generating new release rates is Gaussian: A4$$\begin{array}{rl}g\left(\ln {\mathit{q}}_i\right|{\mathbf{\Lambda }}_i,\ln {\mathit{q}}_j)& =g(\ln q_{k-\mathrm{1}}^{\prime },\ln q_k^{\prime }|\ln q_{k-\mathrm{1}},u_{\ln q})\\ & ={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{\sigma }}}{e}^{-\frac{u_{\ln q}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$ where the new boundary here is generated at hour $k-\mathrm{1}$.

The probability of death of one among the $N_{\mathrm{imp}}-\mathrm{1}$ boundaries is A5$$g\left({\mathbf{\Lambda }}_j\right|{\mathbf{\Lambda }}_i)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{imp}}-\mathrm{1}}},$$ which is a uniform choice among $N_{\mathrm{imp}}-\mathrm{1}$ possibilities. On the other hand, the probability $g\left(\ln {\mathit{q}}_j\right|{\mathbf{\Lambda }}_j,{\mathbf{\Lambda }}_i,\ln {\mathit{q}}_i)$ of destroying the release associated with this position is $\mathrm{1}$.

Data availability

Each dataset used in the paper is available via a DOI on a compliant repository or directly in the cited papers. More precisely, the air concentration observation datasets used in the inversions are available in from http://www.radiochem.org/en/j-online152.html (last access: 3 February 2023) and in , , , and . In addition, the deposition observation dataset used in the inversions is available via 10.5281/zenodo.7016491 . The IBISA (inverse Bayesian inference for source assessment) algorithm, which specifically describes the reversible-jump Markov chain Monte Carlo algorithm used in this article, is available from 10.5281/zenodo.7318543 .

Author contributions

JDLB: software, methodology, conceptualization, investigation, writing – original draft, and visualization. MB: methodology, software, conceptualization, writing – review and editing, visualization, and supervision. OS: resources, methodology, conceptualization, writing – review and editing, visualization, and supervision. YR: methodology, software, and writing – review and editing.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

The authors would like to thank the European Centre for Medium-Range Weather Forecasts (ECMWF) for making available the meteorological fields used in this study. They would also like to thank Sophie Ricci, Lionel Soulhac, Paola Cinnella, Didier Lucor, Yann  Richet, and Anne Mathieu for their suggestions and advice on this work. CEREA is a member of the Institut Pierre-Simon Laplace (IPSL).

Review statement

This paper was edited by Dan Lu and reviewed by two anonymous referees.