QSER Decomposition for Forward and Inverse Modeling of Radon Transport in Lithosphere and Atmosphere

1. Introduction

Radon-222 (${­}^{222}$Rn) is a naturally occurring radioactive gas with a half-life of 3.82 days, continuously produced in soil and rock by the decay of radium-226. Its chemical inertness and short half-life make it an ideal tracer for studying dynamic processes in the lithosphere, atmosphere, and hydrosphere. Applications include earthquake precursor studies [1,2,3], uranium exploration [4], atmospheric transport [5], groundwater-surface water interaction [6], and indoor air quality assessment [7]. The common thread across these applications is the need to interpret measured radon concentrations in terms of physical processes occurring at depth or in the environment.

The mathematical description of radon migration has progressed through several distinct phases. Early models treated radon transport as steady-state Fickian diffusion through homogeneous porous media [8]. These were extended to include advection, radioactive decay, and time-dependent source terms, yielding the standard advection-diffusion-decay equation that remains the foundation of most radon transport models today [4,7]. Numerical implementations enabled simulations in heterogeneous media and complex geometries. Complementing these physical models, statistical and machine learning approaches have been applied to radon time series analysis. Fractal and multifractal methods have revealed long-range correlations not captured by conventional models [1]. Studies have explored the relationship between radon anomalies and seismic activity, demonstrating the utility of radon monitoring for earthquake precursor studies [2]. More recently, artificial neural networks, long short-term memory networks, and support vector regression have been applied to capture nonlinear relationships and improve prediction accuracy. These methods are flexible and data-driven, but they lack physical interpretability.

Despite these advances, several fundamental challenges remain. Models solve directly for the observed concentration, conflating source dynamics with environmental modifications. When a change in radon concentration is detected, it is often unclear whether it originates from the source or from the environment. The medium actively filters the radon response, attenuating amplitude, introducing phase lags, and storing the history of past interactions through diffusive and decay processes. Yet the capacity of the medium to retain and release historical information—its memory—is not treated explicitly. Inverse modeling remains challenging due to ill-posedness and computational cost. Models are fragmented across domains, with different frameworks for soil-to-air, atmospheric, and water transport. Most models lack built-in self-validation, and data integration remains a challenge.

The Source-Environment-Response (SER) framework, introduced by Muhammad et al. [9], addresses several of these challenges by decomposing the radon concentration into a source-driven component and an environmental correction. The SER framework solved the boundary value problem using recorded time-series data as the lower boundary condition, and parameterized the source term through physical quantities [4]. SER advanced the earlier philosophy for solving soil-to-air transport and ion transport from soil to air. However, the SER framework treats the environment as a correction term rather than as a memory accumulator. The memory kernel that governs how past source dynamics influence the current environment remains implicit.

In this paper, we introduce QSER (Generalized Source-Environment-Response), which extends SER by making the memory kernel explicit. The framework is specifically designed for inverse modeling, enabling the recovery of source fields and model parameters from measured concentrations. The paper is organized as follows. Section 2 introduces the QSER framework, followed by it’s application to radon transport in section 3. Section 4 validates QSER against established models from literature, while section 5 explores inverse modeling application via synthetic experiments. Section 6 demonstrates QSER application to inverse modeling in Lithosphere-Atmosphere coupling case, with validation against Nazaroff solution. Section 7 discusses the limitations of QSER and future directions.

2. The QSER Framework

We propose the QSER (Generalized Source-Environment-Response) framework to provide a systematic method for decomposing the observed response of a given dynamical system into two causally distinct components: a conservative source field and an environmental memory field. The framework is built upon an operator split that separates time-reversible (conservative) dynamics from time-irreversible (dissipative) dynamics. The symbol “Q” represents memory of in the QSER framework, this remains consistent throughout the manuscript.

2.1. Physical Motivation for Radon Transport

${­}^{222}Rn$ transport through soil and the atmosphere involves two fundamentally different types of physical processes. Conservative processes, such as advection, transport the radon response without altering its shape or amplitude. These processes are time-reversible and preserve the memory of initial conditions and source history. Dissipative processes, such as molecular diffusion and radioactive decay, smooth gradients, attenuate responses, and remove mass. These processes are time-irreversible and introduce a causal arrow: the present state depends on the entire history of past interactions.

In traditional radon transport models, these two types of processes are coupled in a single equation. The observed concentration $C(\mathbf{r},t)$ is the solution to [9]:

$${\displaystyle \frac{\partial C}{\partial t}}=D{\nabla }^{2}C-\nabla ·\left(\mathbf{v}C\right)-\lambda C+F(\mathbf{r},t)$$

(1)

This equation mixes conservative advection with dissipative diffusion and decay. The QSER framework separates these processes, revealing the distinct roles of source dynamics and environmental memory in shaping the observed radon response.

2.2. Core Decomposition

We define the QSER framework as follows; Consider a system governed by the operator equation:

$$LR=F$$

(2)

where R is the observed response, F is the external forcing or source term, and L is the linear operator describing the system dynamics. The operator L is split into conservative and dissipative parts:

$$L=L_c+L_d$$

(3)

where $L_c$ contains all conservative (time-reversible) terms and $L_d$ contains all dissipative (time-irreversible) terms.

The source field S is defined as the solution to the conservative equation:

$$L_{c}S=F$$

(4)

with the same initial conditions as the observed response R. Physically, S represents what the system would do in the absence of dissipative influences. It preserves the memory of initial conditions and propagates the input F conservatively. For radon transport, S represents the radon concentration that would exist if there were no diffusion or decay — only advective transport.

The environmental field E is defined as the difference between the source field and the observed response:

$$E=S-R$$

(5)

This immediately gives the fundamental decomposition:

$$R=S-E$$

(6)

The environmental field E represents the accumulated history of dissipative interactions: diffusion, decay, boundary effects, and the history of external forcing. It starts from zero initial conditions and grows over time as the system interacts with its environment.

2.3. The Environmental Accumulator Equation

To derive the evolution equation for E, substitute $R=S-E$ into Eq. (2):

$$L(S-E)=F$$

(7)

Expanding $L=L_c+L_d$:

$$L_{c}S+L_{d}S-L_{c}E-L_{d}E=F$$

(8)

Since $L_{c}S=F$ from Eq. (4):

$$F+L_{d}S-L_{c}E-L_{d}E=F$$

(9)

Cancelling F and rearranging:

$$(L_c+L_d)E=L_{d}S$$

(10)

Thus, the environmental field satisfies:

$$LE=L_{d}S$$

(11)

with zero initial conditions. This is the environmental accumulator equation. It states that the environmental field E satisfies the full dissipative operator L, driven by $L_{d}S$ — the dissipative terms acting on the source field.

2.4. Two Limiting Cases

The QSER framework naturally handles two limiting cases.

2.4.1. Pure Conservative Dynamics ($L_d=0$)

When $L_d=0$, the system is purely conservative. The governing equation is $L_{c}R=F$, and the source field satisfies $L_{c}S=F$. Since $L_c=L$, we have $S=R$, and consequently:

$$E=0,R=S$$

(12)

There is no environmental accumulation. The system preserves perfect memory of its initial conditions and the forcing. In radon transport, this corresponds to pure advection without diffusion or decay. This is only possible when radon decay is not significant in the study.

2.4.2. Pure Dissipative Dynamics ($L_c=0$)

When $L_c=0$, the system is purely dissipative. The governing equation is $L_{d}R=F$. The source field does not exist because $L_{c}S=F$ has no solution:

$$S=0$$

(13)

The fundamental decomposition gives:

$$R=-E$$

(14)

Substituting into the governing equation:

$$L_d(-E)=F$$

(15)

$$L_{d}E=-F$$

(16)

The solution for E is:

$$E\left(t\right)=-{\int }_0^{t}Q(t-\tau )F\left(\tau \right)d\tau$$

(17)

and the observed response is:

$$R\left(t\right)={\int }_0^{t}Q(t-\tau )F\left(\tau \right)d\tau$$

(18)

where Q is the memory kernel (Green’s function) of $L_d$. This shows that even in the absence of a conservative core, the system can exhibit memory through the dissipative operator $L_d$ and its memory kernel Q. The response is the accumulated history of the forcing weighted by the memory kernel (Green’s function). In radon transport, this corresponds to pure diffusion and decay without advection — the entire response is environmental.

2.5. Memory Kernel Q (Green’s Function)

For linear systems, the environmental field admits an exact convolution representation using the memory kernel Q. Let $Q(t,\tau )$ satisfy:

$$LQ(t,\tau )=\delta (t-\tau ),Q(t,\tau )=0\mathrm{for}t<\tau$$

(19)

Then the environmental field is given by:

$$E\left(t\right)={\int }_0^{t}Q(t,\tau )L_{d}S\left(\tau \right)d\tau$$

(20)

This representation reveals the environment as a causal memory accumulator. The term $L_{d}S\left(\tau \right)$ is the instantaneous dissipative forcing at time $\tau$, and $Q(t,\tau )$ is the memory kernel that weights how past interactions influence the present environment.

For radon transport in homogeneous media, the memory kernel (Green’s function) takes an exponential form, characteristic of local, short-range memory. The generalization to power-law memory, appropriate for heterogeneous or fractured media, is discussed in Section 8.

2.6. Energy Continuity

In this framework, environment field E interacts with the source field S and dictates how the response is produced. It is therefore, necessary to analyze the nature of this interaction. We in this section discuss the interaction in terms of energy rate of the environment field dynamics. Consider systems where a conserved energy functional exists for $L_c$, we derive a general energy continuity equation for the environmental field. Define the energy functional:

$${\mathcal{E}}_E=\langle E,E\rangle$$

(21)

This definition (21) is well known in fluid flow, Quantum, and signal processing studies. Differentiating with respect to time and using Eq. (11) yields:

$${\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E+{\mathsf{\Phi }}_{\mathrm{boundary}}=J_{S\to E}$$

(22)

where $D_E\ge 0$ is the internal dissipation rate, ${\mathsf{\Phi }}_{\mathrm{boundary}}$ is the energy flux through boundaries, and $J_{S\to E}$ is the energy flux from the source field to the environmental field. Detailed derivation can be seen in the Appendix section. Eq. (22) provides a complete energy accounting: the rate of energy accumulation in E plus dissipation and boundary losses equals the energy transferred from S. Energy is never lost; it is transferred from S to E to irreversible dissipation. A critical insight follows: if $L_d=0$ (no dissipation), then $E\equiv 0$ and no accumulation occurs. Dissipation is not a loss mechanism but the physical process that enables memory accumulation. Without dissipation, the environment remains empty and no causal record is stored.

2.7. Validity Tests

The QSER framework uses four criteria that are system-independent, to ensure and validate the Physics of operator split $L=L_c+L_d$:

2.7.0.1. Test 1: Energy conservation in S.

The source field S must evolve under purely conservative dynamics. In other words, the energy of source dynamics must be conserved:

$${\sigma }_{{\mathcal{E}}_S}<{10}^{-3}$$

(23)

where ${\sigma }_{{\mathcal{E}}_S}$ is the standard deviation of the source energy over the observation window.

2.7.0.2. Test 2: Zero initial condition of E.

The environmental field must start from zero for forward problems. i.e. Problems with known governing equations, seeking to solve for R:

$$E\left(0\right)=0,\dot{E}\left(0\right)=0,\dots$$

(24)

2.7.0.3. Test 3: Exact reconstruction.

The observed response must satisfy $R=S-E$ to within numerical precision:

$$max|R-(S-E)|<{10}^{-3}$$

(25)

2.7.0.4. Test 4: Infinite memory of S (Optional).

The source field must retain perfect memory. Define the unsigned centroid persistence time:

$${\tau }_u={\displaystyle \frac{{\int }_0^{\infty }t|dR/dt|dt}{{\int }_0^{\infty }|dR/dt|dt}}$$

(26)

For a purely conservative system, ${\tau }_u\to \infty$:

$${\tau }_u^S\to \infty \mathrm{as}T\to \infty$$

(27)

The four validity tests are summarized in Table 1. Test 1 ensures the source field remains conservative; Test 2 enforces causality by confirming the environment starts empty; Test 3 verifies the decomposition numerically via exact reconstruction; and Test 4 (optional) confirms the infinite memory of S. Taken together, these provide built-in verification that the operator split has been implemented correctly and that the decomposition preserves the underlying physics. Tests 1–3 are mandatory in any QSER implementation, while Test 4 is recommended when verifying conservative dynamics in well-defined initial-value problems.

3. Radon Transport Equations

The QSER framework is now applied to ${­}^{222}Rn$ transport. We present the general three-dimensional governing equation, perform the operator split, and derive the decomposed system for the source and environmental fields. The solution for linear systems is then presented using the Green’s function (memory kernel) formulation.

3.1. Governing Equation

The transport of ${­}^{222}Rn$ through porous media, fractured rock, and the atmosphere is governed by the advection-diffusion-decay equation. In three-dimensional form, with position vector $\mathbf{r}=(x,y,z)$, the equation is:

$${\displaystyle \frac{\partial C}{\partial t}}=\nabla ·(D\nabla C)-\nabla ·\left(\mathbf{v}C\right)-\lambda C+F(\mathbf{r},t)$$

(28)

where $C(\mathbf{r},t)$ is the radon concentration (Bq m${­}^{-3}$), D is the effective diffusion tensor (m${­}^2$ s${­}^{-1}$), $\mathbf{v}(\mathbf{r},t)$ is the advection velocity vector (m s${­}^{-1}$), $\lambda =2.1\times {10}^{-6}$ s${­}^{-1}$ is the decay constant of ${­}^{222}$Rn, and $F(\mathbf{r},t)$ is the source term representing radon production (Bq m${­}^{-3}$ s${­}^{-1}$). For isotropic media, the diffusion tensor reduces to $D=D\mathbf{I}$, and Eq. (28) becomes:

$${\displaystyle \frac{\partial C}{\partial t}}=D{\nabla }^{2}C-\nabla ·\left(\mathbf{v}C\right)-\lambda C+F$$

(29)

This equation captures the essential physics: molecular diffusion, advective transport, and radioactive decay.

3.2. Operator Split in 3D

Following the QSER framework, we express Eq. (29) in operator form:

$$LR=F$$

(30)

where $R=C$ is the observed response and L is the full operator:

$$L={\displaystyle \frac{\partial }{\partial t}}+\nabla ·(\mathbf{v}·)-D{\nabla }^2+\lambda$$

(31)

The operator is split into conservative and dissipative parts:

$$L=L_c+L_d$$

(32)

The conservative operator $L_c$ contains the time derivative and advection terms, which are time-reversible:

$$L_c={\displaystyle \frac{\partial }{\partial t}}+\nabla ·(\mathbf{v}·)$$

(33)

The dissipative operator $L_d$ contains diffusion and decay, which are time-irreversible:

$$L_d=-D{\nabla }^2+\lambda$$

(34)

This split is physically motivated: advection transports the radon response without altering its shape (conservative), while diffusion smooths gradients and decay removes mass (dissipative).

3.3. Recovering the Source-Environment-Response Framework

The Source-Environment-Response (SER) framework, introduced by Muhammad et al. [9], provides a method for solving the radon transport equation in various media. The SER method decomposes the radon concentration into a source-driven component and an environmental correction [9]:

$$C_{Rn}(x,z,t)={\displaystyle \frac{1}{\lambda }}(S(x,z,t)-E(x,z,t))$$

(35)

The source-driven component $S(x,z,t)$ represents the idealized behavior of radon in the absence of environmental limitations, while the environmental correction $E(x,z,t)$ quantifies deviations caused by real-world constraints such as boundary conditions, dissipation, and local depletion effects [9]. The normalization factor $1/\lambda$ ensures dimensional consistency between production and loss processes. The SER framework was validated against lithosphere-atmosphere coupling [10,11,12,13], radon transport in water [14], and trapped radon in pore channels [9], showing strong agreement with established literature.

The QSER framework generalizes and recovers the SER decomposition as a special case. From the QSER fundamental relation $C=S_{\mathrm{QSER}}-E_{\mathrm{QSER}}$, if we define:

$$S_{\mathrm{QSER}}\equiv {\displaystyle \frac{1}{\lambda }}{S}_{\mathrm{SER}},E_{\mathrm{QSER}}\equiv {\displaystyle \frac{1}{\lambda }}{E}_{\mathrm{SER}}$$

(36)

then:

$$C=S_{\mathrm{QSER}}-E_{\mathrm{QSER}}={\displaystyle \frac{1}{\lambda }}(S_{\mathrm{SER}}-E_{\mathrm{SER}})$$

(37)

Thus, QSER exactly recovers the SER decomposition, with the normalization absorbed into the definitions of the source and environmental fields. While SER provides the core decomposition, QSER extends it by making the memory kernel explicit. In SER, the environmental correction E is defined implicitly; in QSER, E satisfies the environmental accumulator equation $LE=L_{d}S$ and admits the convolution representation $E\left(t\right)={\int }_0^{t}Q(t-\tau )L_{d}S\left(\tau \right)d\tau$, where Q is the memory kernel. This reveals that the environment is not a passive correction term but an active memory accumulator that stores the causal history of dissipative interactions. Additionally, QSER provides four self-validating tests: energy conservation in S, zero initial condition of E, exact reconstruction, and infinite memory of S, which are absent in the original SER framework.

Table 2 summarizes the comparison between SER and QSER.

4. Atmospheric Transport

We now apply the QSER framework to atmospheric radon transport. The full three-dimensional time-dependent solution for the source field is derived from first principles. We begin with the most general form of the atmospheric transport equation and systematically derive the source field S, which we interpret as the ghost of the system — a conservative, memory-preserving component that advects without diffusion or decay. We then generalize the source field to include stochastic velocity and shape-changing behavior, while proving that it remains strictly conservative under these generalizations.

4.1. The General 3D Atmospheric Transport Equation

The transport of a contaminant (including ${­}^{222}Rn$) in the atmosphere is governed by the time-dependent advection-diffusion-decay equation. For a point source release at position ${\mathbf{r}}_0=(0,0,H)$, with general velocity field $\mathbf{v}=(v_x,v_y,v_z)$ and isotropic eddy diffusivity K, the governing equation is:

$${\displaystyle \frac{\partial C}{\partial t}}+\mathbf{v}·\nabla C=K{\nabla }^{2}C-\lambda C+F_0\delta (\mathbf{r}-{\mathbf{r}}_0)\delta \left(t\right)$$

(38)

with zero initial condition:

$$C(\mathbf{r},0)=0$$

(39)

where $C(\mathbf{r},t)$ is the radon concentration (Bq m${­}^{-3}$); $\mathbf{r}=(x,y,z)$ is the position vector; $\mathbf{v}=(v_x,v_y,v_z)$ is the advection velocity vector (m s${­}^{-1}$); K is the effective eddy diffusivity (m${­}^2$ s${­}^{-1}$); $\lambda =2.1\times {10}^{-6}$ s${­}^{-1}$ is the decay constant of ${­}^{222}$Rn; $F_0$ is the total activity released (Bq); ${\mathbf{r}}_0=(0,0,H)$ is the source location; and $\delta (·)$ is the Dirac delta function.

For the special case of wind in the x-direction only, $\mathbf{v}=(u,0,0)$, Eq. (38) reduces to the standard atmospheric dispersion equation:

$${\displaystyle \frac{\partial C}{\partial t}}+u{\displaystyle \frac{\partial C}{\partial x}}=K_y{\displaystyle \frac{{\partial }^{2}C}{\partial y^2}}+K_z{\displaystyle \frac{{\partial }^{2}C}{\partial z^2}}-\lambda C+F_0\delta \left(x\right)\delta \left(y\right)\delta (z-H)\delta \left(t\right)$$

(40)

4.2. Operator Form and Split

Following the QSER framework, we express Eq. (38) in operator form:

$$LC=F$$

(41)

where $R=C$ is the observed response and L is the full operator:

$$L={\displaystyle \frac{\partial }{\partial t}}+\mathbf{v}·\nabla -K{\nabla }^2+\lambda$$

(42)

and the forcing term is:

$$F=F_0\delta (\mathbf{r}-{\mathbf{r}}_0)\delta \left(t\right)$$

(43)

The operator L is split into conservative and dissipative parts:

$$L=L_c+L_d$$

(44)

The conservative operator $L_c$ contains the time derivative and advection terms, which are time-reversible:

$$L_c={\displaystyle \frac{\partial }{\partial t}}+\mathbf{v}·\nabla$$

(45)

The dissipative operator $L_d$ contains diffusion and decay, which are time-irreversible:

$$L_d=-K{\nabla }^2+\lambda$$

(46)

This split is physically motivated: advection transports the radon response without altering its shape (conservative), while diffusion smooths gradients and decay removes mass (dissipative).

4.3. The Source Field S — The Ghost

The source field S is defined as the solution to the conservative equation:

$$L_{c}S=F$$

(47)

with the same initial condition as the observed concentration:

$$S(\mathbf{r},0)=0$$

(48)

Substituting Eqs. (45) and (43):

$${\displaystyle \frac{\partial S}{\partial t}}+\mathbf{v}·\nabla S=F_0\delta (\mathbf{r}-{\mathbf{r}}_0)\delta \left(t\right)$$

(49)

$$S(\xi ,t)=F_{0}H\left(t\right)\delta (\xi -{\mathbf{r}}_0)$$

(50)

where $H\left(t\right)$ is the Heaviside step function enforcing causality, and $\xi =\mathbf{r}-\mathbf{v}t$:

$$S(\mathbf{r},t)=F_{0}H\left(t\right)\delta (\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t)$$

(51)

4.3.1. Physical Interpretation

The source field S in Eq. (51) represents the radon concentration that would exist in the absence of diffusion and decay. It has the following properties:

1.

Pure advection: The delta function moves with velocity $\mathbf{v}$.

2.

No diffusion: The delta function remains sharp (no spreading).

3.

No decay: There is no exponential attenuation factor $e^{-\lambda t}$.

4.

Perfect memory: The source retains the history of its initial condition.

5.

Causality: The Heaviside function $H\left(t\right)$ ensures no response before $t=0$.

4.4. The Generalized Source Field — Stochastic and Shape-Changing

We now generalize the source field to include:

1.

Stochastic velocity (random path),

2.

Shape-changing behavior (squeezing, expanding, stretching),

3.

A general normalized shape function.

4.4.1. General Form

The generalized source field is written as:

$$S(\mathbf{r},t)=F_{0}H\left(t\right)f\left(\mathbf{r}-\xi \left(t\right),t\right)$$

(52)

where $\xi \left(t\right)$ is the trajectory of the ghost; and $f(\eta ,t)$ is the shape function.

4.4.2. The Trajectory

The trajectory $\xi \left(t\right)$ is defined by:

$$\xi \left(t\right)={\mathbf{r}}_0+{\int }_0^t\mathbf{v}\left(\tau \right)d\tau$$

(53)

where $\mathbf{v}\left(\tau \right)$ can be: deterministic, $\mathbf{v}\left(\tau \right)=\overline{\mathbf{v}}$; stochastic, $\mathbf{v}\left(\tau \right)=\overline{\mathbf{v}}+\sigma {\xi }_v\left(\tau \right)$; or any integrable function (deterministic or stochastic).

4.4.3. The Shape Function

The shape function $f(\eta ,t)$ must satisfy the normalization condition:

$${\int }_{{\mathbb{R}}^3}f(\eta ,t)d\eta =1\mathrm{for}\mathrm{all}t$$

(54)

This condition ensures that the ghost remains conservative regardless of shape changes.

4.4.4. Examples of Shape Functions

4.4.4.1. Gaussian (Expanding/Squeezing):

$$f(\eta ,t)={\displaystyle \frac{1}{{\left(2\pi \sigma {\left(t\right)}^2\right)}^{3/2}}}exp\left(-{\displaystyle \frac{{\left|\eta \right|}^2}{2\sigma {\left(t\right)}^2}}\right)$$

(55)

where $\sigma \left(t\right)$ controls the width: constant $\sigma$ yields a fixed width; increasing $\sigma \left(t\right)$ yields an expanding ghost; decreasing $\sigma \left(t\right)$ yields a squeezing ghost; and oscillating $\sigma \left(t\right)$ yields a breathing ghost.

4.4.4.2. Rectangular (Expanding/Squeezing):

$$f(\eta ,t)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{V\left(t\right)}},\hfill & \left|\eta \right|<a\left(t\right)\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.,V\left(t\right)={\displaystyle \frac{4}{3}}\pi a{\left(t\right)}^3$$

(56)

4.4.4.3. General Normalized Function:

Any function satisfying Eq. (54) is admissible.

4.4.5. Conservative Property

We now prove that the generalized source field S in Eq. (52) is conservative.

Proof. 

The total mass of the ghost is:

$$M_S\left(t\right)={\int }_{{\mathbb{R}}^3}S(\mathbf{r},t)d\mathbf{r}$$

(57)

Substituting Eq. (52):

$$M_S\left(t\right)=F_{0}H\left(t\right){\int }_{{\mathbb{R}}^3}f(\mathbf{r}-\xi \left(t\right),t)d\mathbf{r}$$

(58)

Let $\eta =\mathbf{r}-\xi \left(t\right)$. Then $d\mathbf{r}=d\eta$:

$$M_S\left(t\right)=F_{0}H\left(t\right){\int }_{{\mathbb{R}}^3}f(\eta ,t)d\eta$$

(59)

Using the normalization condition in Eq. (54):

$$M_S\left(t\right)=F_{0}H\left(t\right)$$

(60)

Therefore:

$${\displaystyle \frac{d{M}_S}{dt}}=0\mathrm{for}t>0$$

(61)

The total mass of the ghost is conserved regardless of:

• Whether $\mathbf{v}\left(t\right)$ is deterministic or stochastic,
• Whether $f(\eta ,t)$ changes shape.

□

4.4.6. Special Cases

Table 3 summarizes the special cases of the generalized source field.

4.5. Numerical Representation of the Source Field

In numerical implementations, the Dirac delta function in Eq. (51) cannot be represented directly on a discrete grid. We therefore approximate it using a Gaussian:

$$\delta (\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t)\approx {\displaystyle \frac{1}{{\left(2\pi {\sigma }^2\right)}^{3/2}}}exp\left(-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}_0{-\mathbf{v}t|}^2}{2{\sigma }^2}}\right)$$

(62)

where $\sigma$ is a small parameter controlling the width of the numerical delta. The approximation becomes exact in the limit $\sigma \to 0$.

The Gaussian approximation for the source field is therefore:

$$S(\mathbf{r},t)\approx {\displaystyle \frac{F_{0}H\left(t\right)}{{\left(2\pi {\sigma }^2\right)}^{3/2}}}exp\left(-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}_0{-\mathbf{v}t|}^2}{2{\sigma }^2}}\right)$$

(63)

For the generalized stochastic shape-changing source field:

$$S(\mathbf{r},t)\approx {\displaystyle \frac{F_{0}H\left(t\right)}{{\left(2\pi \sigma {\left(t\right)}^2\right)}^{3/2}}}exp\left(-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}_0-{\int }_0^t{\mathbf{v}\left(\tau \right)d\tau |}^2}{2\sigma {\left(t\right)}^2}}\right)$$

(64)

4.6. Summary and Transition

In this section, we have fully characterized the source field S:

$$S(\mathbf{r},t)=F_{0}H\left(t\right)f\left(\mathbf{r}-{\mathbf{r}}_0-{\int }_0^t\mathbf{v}\left(\tau \right)d\tau ,t\right)$$

(65)

with the conservative condition:

$${\int }_{{\mathbb{R}}^3}f(\eta ,t)d\eta =1\mathrm{for}\mathrm{all}t$$

(66)

Key properties established:

1.

The source field is strictly conservative.

2.

It admits stochastic velocity without losing conservativity.

3.

It admits shape-changing behavior without losing conservativity.

4.

It represents the “ghost” of the system — invisible yet ever-present.

Figure 1 illustrates four realizations of the generalized source field at $t=3.0$ s. The total mass of the source field is $M_S\left(t\right)={\int }_{{\mathbb{R}}^3}S(\mathbf{r},t)d\mathbf{r}=F_{0}H\left(t\right)$, which follows directly from the normalization condition. Therefore $d{M}_S/dt=0$ for $t>0$, proving that the source field remains strictly conservative regardless of whether the velocity is deterministic or stochastic and regardless of whether the shape function changes over time. This conservative property has significant implications for inverse modeling. When observed concentration data deviate from the source field behavior, the difference is attributed entirely to the environmental field E, which accumulates the history of dissipative interactions. The generalized source field therefore provides a flexible yet physically constrained baseline for interpreting radon transport data, where velocity fields may be uncertain or source shapes may evolve due to changing environmental conditions. The ability to accommodate stochasticity and shape-changing behavior without compromising mass conservation makes the QSER framework suitable for real-world applications where these complexities are present.

5. Memory Kernel and Environmental Field

Having fully characterized the source field S in previous, we now derive the memory kernel Q and the environmental field E. Together with S, these complete the QSER decomposition and yield the full solution for the observed concentration C.

5.1. The Memory Kernel Q

The memory kernel Q is the Green’s function of the full operator L. It satisfies:

$$LQ(\mathbf{r},t;{\mathbf{r}}^{\prime },\tau )=\delta (\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-\tau )$$

(67)

with:

• Causality: $Q=0$ for $t<\tau$,
• Boundary condition: $Q\to 0$ as $\left|\mathbf{r}\right|\to \infty$.

$$Q(\mathbf{r},t;{\mathbf{r}}^{\prime },\tau )={\displaystyle \frac{H(t-\tau )}{{\left(4\pi K(t-\tau )\right)}^{3/2}}}exp\left(-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}^{\prime }{-\mathbf{v}(t-\tau )|}^2}{4K(t-\tau )}}-\lambda (t-\tau )\right)$$

(68)

The memory kernel $Q\left(s\right)$ governs how past emissions influence the present concentration. For a source emission at time $\tau$, the influence at time t decays as the time lag $s=t-\tau$ increases. Figure 2 illustrates this decay for different diffusion coefficients K.

The left panel shows $Q\left(s\right)$ at the source position. The kernel decays rapidly for larger K because diffusion spreads the memory over a larger spatial extent, reducing the peak amplitude. The decay is governed by both the diffusive spreading term ${\left(4\pi Ks\right)}^{-3/2}$ and the radioactive decay term $e^{-\lambda s}$. For ${­}^{222}Rn$, $\lambda =2.1\times {10}^{-6}$ s${­}^{-1}$ is small, so the decay is dominated by diffusion over the time scales shown.

The right panel shows $Q\left(s\right)$ at a fixed detector position $(5,3)$ m. The detector does not see the emission immediately. Instead, the kernel shows a peak at the advective travel time, with the arrival time and peak amplitude depending on K. For smaller K, the signal arrives later and is more concentrated. For larger K, diffusion spreads the signal, causing earlier arrival and reduced peak amplitude. The detector response represents what a measurement device would record from a past emission, with the convolution over all past emissions yielding the observed concentration.

The memory kernel therefore determines how the system forgets. It writes the history of past emissions into the environmental field E, which accumulates the integrated history of Q weighted by the dissipative source term $L_{d}S$. In the discussion, we note that $\lambda$ may represent reaction rates, chemical loss, or other decay processes in applications beyond radon transport.

5.2. The Environmental Field E

The environmental field E satisfies the accumulator equation:

$$LE=L_{d}S$$

(69)

with zero initial condition:

$$E(\mathbf{r},0)=0$$

(70)

5.2.1. Solution Using Convolution

Since Q is the Green’s function of L, the solution for E is:

$$E(\mathbf{r},t)={\int }_0^t{\int }_{{\mathbb{R}}^3}Q(\mathbf{r},t;{\mathbf{r}}^{\prime },\tau )L_{d}S({\mathbf{r}}^{\prime },\tau )d{\mathbf{r}}^{\prime }d\tau$$

(71)

Substituting the explicit form of Q from Eq. (68), we obtain:

$$E(\mathbf{r},t)=F_0{\int }_0^t{\displaystyle \frac{1}{{\left(4\pi Ks\right)}^{3/2}}}exp\left(-{\displaystyle \frac{{\left|\xi \right|}^2}{4Ks}}-\lambda s\right)\left[-{\displaystyle \frac{{\left|\xi \right|}^2}{4K{s}^2}}+{\displaystyle \frac{3}{2s}}+\lambda \right]ds$$

(72)

where:

$$\xi =\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t$$

(73)

and $s=t-\tau$.

Figure 3 illustrates the environmental field for different diffusion coefficients. The left panel shows $E\left(s\right)$ at the source position over the early time lag range where accumulation is most dynamic. The environment grows from zero, initially dominated by diffusion effects (positive E), before transitioning through zero as decay effects become significant (negative E). The diffusion coefficient K controls both the magnitude and the rate of this transition. Larger K produces faster accumulation and earlier crossover to negative values, while smaller K produces slower accumulation and delayed crossover. The right panel shows a 2D visualization of E at $t=3.0$ s with the source trajectory overlaid. The environmental field forms a shadow along the source path, representing the accumulated history of dissipative interactions. The observed concentration C is recovered through the fundamental decomposition $C=S-E$, meaning what is measured is the conservative source minus the accumulated environmental memory.

5.3. The Full QSER Solution

The observed concentration C is recovered through the fundamental QSER decomposition:

$$C(\mathbf{r},t)=S(\mathbf{r},t)-E(\mathbf{r},t)$$

(74)

5.3.1. Complete Solution

Substituting Eqs. (51) and (72):

$$C(\mathbf{r},t)=F_{0}H\left(t\right)\delta (\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t)-F_0{\int }_0^t{\displaystyle \frac{1}{{\left(4\pi Ks\right)}^{3/2}}}exp\left(-{\displaystyle \frac{{\left|\xi \right|}^2}{4Ks}}-\lambda s\right)\left[-{\displaystyle \frac{{\left|\xi \right|}^2}{4K{s}^2}}+{\displaystyle \frac{3}{2s}}+\lambda \right]ds$$

(75)

where:

$$\xi =\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t$$

(76)

5.4. Special Cases and Reductions

5.4.1. Case 1: Pure Advection ($K=0,\lambda =0$)

When diffusion and decay are absent:

$$C(\mathbf{r},t)=F_{0}H\left(t\right)\delta (\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t)$$

(77)

The environment is zero, and the ghost is visible.

5.4.2. Case 2: Pure Diffusion ($\mathbf{v}=0,\lambda =0$)

When advection and decay are absent:

$$C(\mathbf{r},t)=-F_0{\int }_0^t{\displaystyle \frac{1}{{\left(4\pi Ks\right)}^{3/2}}}exp\left(-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}_0{|}^2}{4Ks}}\right)\left[-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}_0{|}^2}{4K{s}^2}}+{\displaystyle \frac{3}{2s}}\right]ds$$

(78)

5.4.3. Case 3: Steady-State Limit ($t\to \infty$)

In the limit $t\to \infty$, the source term vanishes and the environmental field approaches the steady-state Gaussian plume. Evaluating the integral in Eq. (72):

$$C\left(\mathbf{r}\right)={\displaystyle \frac{F_0}{4\pi |\mathbf{r}-{\mathbf{r}}_0|K}}exp\left(-{\displaystyle \frac{|\mathbf{r}-{\mathbf{r}}_0{|}^2}{4K|\mathbf{r}-{\mathbf{r}}_0|/|\mathbf{v}|}}-{\displaystyle \frac{\lambda |\mathbf{r}-{\mathbf{r}}_0|}{\left|\mathbf{v}\right|}}\right)$$

(79)

5.4.4. Case 4: Standard Gaussian Plume

For ground-level concentration ($z=0$) with wind in the x-direction only ($\mathbf{v}=(u,0,0)$), and using the dispersion coefficients:

$${\sigma }_y^2={\displaystyle \frac{2K_{y}x}{u}},{\sigma }_z^2={\displaystyle \frac{2K_{z}x}{u}}$$

(80)

Eq. (79) reduces to the classical Gaussian plume:

$$C(x,y,0)={\displaystyle \frac{F_0}{\pi u{\sigma }_y{\sigma }_z}}exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}-{\displaystyle \frac{\lambda x}{u}}\right)$$

(81)

5.5. Recovery of Industry-Standard Atmospheric Dispersion Models

The QSER framework recovers all major industry-standard atmospheric dispersion models as special cases. These models fall into two categories: straight-line Gaussian plumes and sector-averaged Gaussian plumes.

5.5.1. Straight-Line Gaussian Plume Models

For models that do not average over wind sectors, the QSER steady-state solution reduces directly to the classical Gaussian plume. For a source at height H with wind in the x-direction only, the ground-level concentration at $z=0$ with reflection is:

$$C(x,y,0)={\displaystyle \frac{F_0}{\pi u{\sigma }_y{\sigma }_z}}exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}\right)exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)$$

(82)

This is exactly the form used in ALOHA and CAMEO emergency response models [18], as well as standard Gaussian plume formulations [7]. The dispersion coefficients ${\sigma }_y$ and ${\sigma }_z$ are determined by the atmospheric stability class and downwind distance x, as parameterized in each model.

5.5.2. Sector-Averaged Gaussian Plume Models

Models such as RESRAD-OFFSITE [15], CAP88-PC [16], ISCLT3, and SCREEN3 [17] average the concentration over a $22.5^{\circ }$ wind sector. The sector width at downwind distance x is $W=2xtan(11.{25}^{\circ })$. The sector-averaged concentration is:

$$\overline{C}\left(x\right)={\displaystyle \frac{1}{W}}{\int }_{-W/2}^{W/2}C(x,y,0)dy$$

(83)

Substituting Eq. (82) and integrating over the sector width (where the sector width is large compared to the plume width, allowing the integral to extend to $\pm \infty$) gives:

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0}{\sqrt{2\pi }u{\sigma }_z\left(x\right)xtan(11.{25}^{\circ })}}exp\left[-{\displaystyle \frac{H^2}{2{\sigma }_z^2\left(x\right)}}-{\displaystyle \frac{\lambda x}{u}}\right]$$

(84)

This exactly matches the RESRAD-OFFSITE sector-averaged Gaussian plume equation [15]. The full RESRAD-OFFSITE formulation includes source depletion from wet and dry deposition:

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0}{\sqrt{2\pi }{U}_{\mathrm{eff}}{\sigma }_z\left(x\right)xtan(11.{25}^{\circ })}}exp\left[-{\displaystyle \frac{H_{\mathrm{eff}}^2}{2{\sigma }_z^2\left(x\right)}}\right]exp\left[-{\int }_0^{x}F\left(s\right)ds\right]$$

(85)

where, $U_{\mathrm{eff}}$ is the effective wind speed, $H_{\mathrm{eff}}$ is the effective release height, and $F\left(s\right)$ is the depletion factor. CAP88-PC [16] uses the same sector-averaged formulation with differences only in the parameterization of vertical dispersion coefficients and plume rise. ISCLT3 and SCREEN3 [17] are also recovered by reduction to the same Gaussian plume foundation. Table 4 summarizes the industry-standard atmospheric models recovered by QSER. The full derivation is provided in Appendix A.

5.5.3. What QSER Adds Beyond Industry Models

While QSER recovers these industry models as special cases, it provides several additional capabilities not present in the standard formulations. The decomposition $C=S-E$ explicitly separates what arrives from the source from what the environment has done to it, providing physical interpretability. The memory kernel Q weights the historical influence of past emissions on current concentrations, enabling time-dependent analysis that is absent in steady-state models. The energy continuity equation ${\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E=J_{S\to E}$ provides quantitative accounting of energy flow from source to environment to dissipation. QSER provides the full time-dependent solution, whereas most industry models are steady-state. Finally, QSER allows the source field S to be solved analytically while the environmental field E can be solved numerically, enabling flexible implementation for inverse modeling applications where source identification and parameter estimation are required [19].

6. Lithosphere-Atmosphere Transport

To demonstrate the applicability of QSER in lithosphere-atmosphere transport, we compare with the classical Nazaroff solution [7] for steady-state radon transport in homogeneous soil:

$$D{\displaystyle \frac{d^{2}C}{d{z}^2}}-v{\displaystyle \frac{dC}{dz}}-\lambda C+F=0$$

(86)

with boundary conditions $C\left(0\right)=0$ and $C(\infty )=C_{\infty }$. The Nazaroff solution is:

$$C_{\mathrm{Nazaroff}}\left(z\right)=C_{\infty }\left(1-e^{-az}\right),C_{\infty }={\displaystyle \frac{F}{\lambda }},a={\displaystyle \frac{v}{2D}}+\sqrt{{\displaystyle \frac{v^2}{4D^2}}+{\displaystyle \frac{\lambda }{D}}}$$

(87)

Applying the QSER decomposition $C=S-E$, the source field satisfies the conservative advection equation $-vdS/dz=F$ with $S\left(0\right)=0$, yielding:

$$S\left(z\right)={\displaystyle \frac{F}{v}}z$$

(88)

The environmental field then follows as $E=S-C_{\mathrm{Nazaroff}}$:

$$E\left(z\right)={\displaystyle \frac{F}{v}}z-{\displaystyle \frac{F}{\lambda }}\left(1-e^{-az}\right)$$

(89)

The QSER decomposition enables complete recovery of all physical parameters from the decomposed fields: F from $C_{\infty }$, v from the slope of S, a from the exponential tail of E, and D from the relationship $D=v/a-\lambda /a^2$.

Figure 4 presents the QSER decomposition and parameter recovery. Panel (a) confirms that $C=S-E$ exactly reproduces the Nazaroff profile, validating the decomposition. Panel (b) reveals the physical separation: the ghost field $S\left(z\right)$ (red) is strictly linear, representing purely advective transport without dissipation — it carries the memory of the source forward. The environmental field $E\left(z\right)$ (green) starts from zero at the surface and accumulates the history of diffusive and decay effects. Panel (c) confirms the linearity of S with slope $F/v=100,000$, demonstrating its conservative nature. Panel (d) shows the exponential component $E_{\exp }\left(z\right)=C_{\infty }{e}^{-az}$ (blue points) with exponential fit (red dashed) yielding $a=0.2002$ m${­}^{-1}$.

The parameter recovery is summarized in Table 5.

All parameters are recovered with near-perfect accuracy (errors $<0.21\%$), confirming that QSER preserves the complete physics of the transport problem. The 0.21% error in D is within numerical precision and arises from the propagation of floating-point errors in the recovery chain.

The QSER framework provides several insights beyond the traditional Nazaroff solution. First, it reveals that the observed concentration is not a fundamental quantity but the difference between two competing fields — the ghost (conservative memory) and the environment (dissipative accumulator). Second, it shows that the environment actively stores the history of past source dynamics through the memory kernel Q, making explicit what traditional models treat implicitly. Third, the recovery of all physical parameters from the decomposed fields demonstrates that QSER has powerful inverse modeling capacity: given observed concentration data, one can estimate the source parameters by finding the ghost field S that best satisfies $S=C_{\mathrm{obs}}+E$, with E governed by the environmental accumulator equation $LE=L_{d}S$. This transforms the ill-posed inverse problem into a well-posed optimization over conservative source parameters only.

The same decomposition applies to atmospheric plumes, where the ghost field $S(\mathbf{r},t)=F_{0}H\left(t\right)\delta (\mathbf{r}-{\mathbf{r}}_0-\mathbf{v}t)$ represents purely advective transport, and the environmental field accumulates the history of diffusion and decay via the memory kernel Q. In the steady-state limit, this recovers the classical Gaussian plume, demonstrating that QSER unifies lithosphere and atmosphere transport under a single rigorous framework.

7. Remarks

The QSER framework introduced in this paper provides a new perspective for studying radon transport across lithosphere and atmosphere domains. By decomposing the observed concentration into a conservative ghost field S and an environmental memory field E, QSER reveals the fundamental duality in radon transport: the ghost carries the memory of the source forward without loss, while the environment accumulates the history of dissipative interactions through the memory kernel Q.

The key contributions of this work are:

1.

Explicit memory kernel: Unlike the SER framework [9] which treated the environment as a passive correction, QSER makes the memory kernel Q explicit through the environmental accumulator equation $LE=L_{d}S$. This reveals that the environment is not a loss mechanism but an active memory accumulator that stores causal history.

2.

Operator splitting: The separation of the transport operator into conservative ($L_c$) and dissipative ($L_d$) parts provides physical interpretability. The ghost field S satisfies purely conservative dynamics while E captures all dissipative effects.

3.

Generalized source field: The source field is extended to include stochastic velocity fields and shape-changing behavior while preserving strict mass conservation through the normalization condition $\int f(\eta ,t)d\eta =1$.

4.

Self-validation: Four system-independent tests (energy conservation in S, zero initial condition of E, exact reconstruction $C=S-E$, and infinite memory of S) provide built-in verification of the decomposition.

5.

Recovery of industry models: QSER recovers all major atmospheric dispersion models (Gaussian plume, RESRAD-OFFSITE, CAP88-PC, ISCLT3, SCREEN3, ALOHA/CAMEO) as special cases while providing additional time-dependent capabilities.

6.

Lithosphere-atmosphere unification: The same QSER decomposition applies to both soil radon profiles (Nazaroff solution) and atmospheric plumes, providing a consistent language for radon transport across domains.

The validation against the Nazaroff solution demonstrates near-perfect recovery of all physical parameters (F, v, a, D) with errors below $0.21\%$. This confirms that QSER preserves the complete physics of the transport problem. The parameter recovery framework transforms the ill-posed inverse problem into a well-posed sequential optimization over conservative source parameters: F from $C_{\infty }$, v from the slope of S, a from the exponential tail of E, and D from the relationship $D=v/a-\lambda /a^2$.

All validations in this paper are performed against established literature models that have been extensively applied to real field data. The Nazaroff solution [7] has been validated against soil radon measurements worldwide. The Gaussian plume and sector-averaged models (RESRAD-OFFSITE [15], CAP88-PC [16], ISCLT3 [17]) are regulatory standards used in real-world environmental monitoring and emergency response. The memory kernel and environmental field behavior have been qualitatively confirmed against atmospheric radon observations [5] and soil radon time series [10,11,12,13]. By extension, the QSER framework inherits the practical applicability of these validated models while providing additional interpretive and analytical capabilities not present in the original formulations.

The ability to separate source dynamics from environmental memory has significant implications for real-world applications. In earthquake precursor studies [2,3], distinguishing between source-driven anomalies (changes in radon production) and environmentally-driven variations (changes in diffusion due to permeability changes) is critical. QSER provides a systematic framework for this separation. In atmospheric monitoring, identifying whether a detected radon peak originates from a new release or from environmental accumulation is essential for source apportionment.

Future work includes: (1) application to real radon time series data with Monte Carlo uncertainty quantification using the Metropolis-Hastings algorithm [20,21], (2) extension to power-law memory for heterogeneous and fractured media where the memory kernel takes the form $Q\left(s\right)\propto s^{-\alpha }$ with $0<\alpha <1$, (3) integration with physics-informed neural networks [19] for inverse modeling where the physics loss is the QSER decomposition itself, (4) development of FQSER (fractional QSER) for anomalous diffusion regimes where the diffusion operator is fractional, and (5) open-source code release to enable community adoption and application to diverse radon datasets.

The QSER framework opens new possibilities for source apportionment, plume fingerprinting, and real-time environmental monitoring where separating source dynamics from environmental memory is essential. By making the memory kernel explicit and providing self-validating tests, QSER addresses the fundamental challenges identified in the introduction: the conflation of source and environment, the absence of explicit memory treatment, the ill-posedness of inverse problems, and the fragmentation of models across domains.

Appendix A.1. The QSER Gaussian Plume Solution

Starting from the QSER steady-state solution for a point source at height H with wind in the x-direction only, the ground-level concentration at $z=0$ with reflection is:

$$C(x,y,0)={\displaystyle \frac{F_0}{\pi u{\sigma }_y{\sigma }_z}}exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}\right)exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)$$

(A1)

where $F_0$ is the source strength, u is the wind speed, ${\sigma }_y$ and ${\sigma }_z$ are the lateral and vertical dispersion coefficients, H is the effective release height, and $\lambda$ is the decay constant.

Appendix A.2. Sector Averaging

For a sector of angular width $22.5^{\circ }$, the sector width at downwind distance x is:

$$W=2xtan(11.{25}^{\circ })$$

(A2)

The sector-averaged concentration is defined as:

$$\overline{C}\left(x\right)={\displaystyle \frac{1}{W}}{\int }_{-W/2}^{W/2}C(x,y,0)dy$$

(A3)

Appendix A.3. Substitution and Integration

Substitute Eq. (A1) into Eq. (A3):

$$\overline{C}\left(x\right)={\displaystyle \frac{1}{W}}{\displaystyle \frac{F_0}{\pi u{\sigma }_y{\sigma }_z}}exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right){\int }_{-W/2}^{W/2}exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}\right)dy$$

(A4)

For sector-averaged models, the sector width W is large compared to the plume width ${\sigma }_y$ [15]. Therefore:

$${\int }_{-W/2}^{W/2}exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}\right)dy\approx {\int }_{-\infty }^{\infty }exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}\right)dy$$

(A5)

The Gaussian integral is:

$${\int }_{-\infty }^{\infty }exp\left(-{\displaystyle \frac{y^2}{2{\sigma }_y^2}}\right)dy=\sqrt{2\pi }{\sigma }_y$$

(A6)

Appendix A.4. Simplification

Substitute Eq. (A6) into Eq. (A4):

$$\overline{C}\left(x\right)={\displaystyle \frac{1}{W}}{\displaystyle \frac{F_0}{\pi u{\sigma }_y{\sigma }_z}}exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)\sqrt{2\pi }{\sigma }_y$$

(A7)

Cancel ${\sigma }_y$:

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0\sqrt{2\pi }}{\pi Wu{\sigma }_z}}exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)$$

(A8)

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0\sqrt{2}}{\sqrt{\pi }Wu{\sigma }_z}}exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)$$

(A9)

Appendix A.5. Substitute Sector Width

Substitute $W=2xtan(11.{25}^{\circ })$:

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0\sqrt{2}}{\sqrt{\pi }2xtan(11.{25}^{\circ })u{\sigma }_z}}exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)$$

(A10)

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0}{\sqrt{2\pi }u{\sigma }_z\left(x\right)xtan(11.{25}^{\circ })}}exp\left(-{\displaystyle \frac{H^2}{2{\sigma }_z^2\left(x\right)}}\right)exp\left(-{\displaystyle \frac{\lambda x}{u}}\right)$$

(A11)

Appendix A.6. Comparison with RESRAD-OFFSITE

RESRAD-OFFSITE [15] gives the sector-averaged concentration as:

$$\overline{C}\left(x\right)={\displaystyle \frac{F_0}{\sqrt{2\pi }{U}_{\mathrm{eff}}{\sigma }_z\left(x\right)xtan(11.{25}^{\circ })}}exp\left(-{\displaystyle \frac{H_{\mathrm{eff}}^2}{2{\sigma }_z^2\left(x\right)}}\right)exp\left(-{\int }_0^{x}F\left(s\right)ds\right)$$

(A12)

where $U_{\mathrm{eff}}$ is the effective wind speed, $H_{\mathrm{eff}}$ is the effective release height, and $F\left(s\right)$ is the depletion factor accounting for wet and dry deposition.

Appendix A.7. Verification

Comparing Eq. (A11) and Eq. (A12), the mathematical structure is identical. The differences are parameterization conventions: u versus $U_{\mathrm{eff}}$, H versus $H_{\mathrm{eff}}$, and $\lambda x/u$ versus $\int F\left(s\right)ds$. The QSER framework therefore exactly recovers the RESRAD-OFFSITE sector-averaged Gaussian plume formulation [15]. CAP88-PC [16] uses the same formulation with different parameterizations of ${\sigma }_z$ and plume rise. ISCLT3 and SCREEN3 [17] are recovered by reduction to the same Gaussian plume foundation. For straight-line models without sector averaging (ALOHA, CAMEO [18], standard Gaussian plume), QSER recovers the direct Gaussian plume solution in Eq. (A1) without the sector integration step.

Appendix B.1. The Operator Formulation

Let the observed response R satisfy the governing equation:

$$LR=F$$

(A13)

where L is the full system operator and F is the external forcing. The operator L is split into conservative and dissipative parts:

$$L=L_c+L_d$$

(A14)

where $L_c$ contains all conservative (time-reversible) terms; and $L_d$ contains all dissipative (time-irreversible) terms.

Appendix B.2. Definition of Source and Environmental Fields

The source field S is defined as the solution to the conservative equation:

$$L_{c}S=F$$

(A15)

with the same initial conditions as the observed response R.

The environmental field E is defined as the difference:

$$E\equiv S-R$$

(A16)

From Eqs. (A13) and (A15), we can derive the governing equation for E:

$$\begin{array}{cc}\hfill L(S-E)& =F\hfill \\ \hfill L_{c}S+L_{d}S-L_{c}E-L_{d}E& =F\hfill \\ \hfill F+L_{d}S-L_{c}E-L_{d}E& =F\hfill \\ \hfill (L_c+L_d)E& =L_{d}S\hfill \end{array}$$

(A17)

Thus:

$$\begin{array}{|c|}\hline LE=L_{d}S\\ \hline\end{array}$$

(A18)

with zero initial conditions: $E\left(0\right)=0$, $\dot{E}\left(0\right)=0$, etc.

Appendix B.3. Derivation of the Energy Continuity Equation

We now derive the energy continuity equation for the environmental field. The derivation proceeds in several steps.

Appendix B.3.1. Step 1: Power Pairing

Take the inner product (power pairing) of Eq. (A18) with $\dot{E}$:

$$\langle \dot{E},LE\rangle =\langle \dot{E},L_{d}S\rangle$$

(A19)

where $\langle ·,·\rangle$ denotes an appropriate inner product or duality pairing. For ODEs, this is simply multiplication; for PDEs, it represents integration over the spatial domain.

Appendix B.3.2. Step 2: Split the Operator

Substitute $L=L_c+L_d$ into Eq. (A19):

$$\langle \dot{E},L_{c}E\rangle +\langle \dot{E},L_{d}E\rangle =\langle \dot{E},L_{d}S\rangle$$

(A20)

Appendix B.3.3. Step 3: Identify the Energy Rate

For a conservative operator $L_c$, the power supplied by the conservative dynamics equals the rate of change of stored energy:

$${\displaystyle \frac{d{\mathcal{E}}_E}{dt}}\equiv \langle \dot{E},L_{c}E\rangle$$

(A21)

This is the definition of the energy functional ${\mathcal{E}}_E$ in the context of the environmental field. It is not an assumption but a statement of what we mean by energy.

Appendix B.3.4. Step 4: Define Dissipation and Flux

Define the dissipation term:

$$D_E\equiv \langle \dot{E},L_{d}E\rangle$$

(A22)

Define the source-to-environment energy flux:

$$J_{S\to E}\equiv \langle \dot{E},L_{d}S\rangle$$

(A23)

Appendix B.3.5. Step 5: The Energy Continuity Equation

Substituting Eqs. (A21), (A22), and (A23) into Eq. (A20) yields the energy continuity equation:

$$\begin{array}{|c|}\hline {\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E=J_{S\to E}\\ \hline\end{array}$$

(A24)

Appendix B.4. Variational Derivation

The energy continuity equation can also be derived from a variational principle. Consider the Lagrangian:

$$\mathcal{L}={\displaystyle \frac{1}{2}}\langle E,LE\rangle -\langle E,L_{d}S\rangle$$

(A25)

The Euler-Lagrange equations yield:

$$LE=L_{d}S$$

(A26)

The energy functional is obtained from the Legendre transform:

$${\mathcal{E}}_E=\langle \pi ,E\rangle -\mathcal{L}$$

(A27)

where $\pi =\partial \mathcal{L}/\partial \dot{E}$ is the canonical momentum. For conservative systems, this yields the standard energy expression.

Appendix B.5. Boundary Terms for PDEs

For systems with spatial boundaries, the power pairing in Eq. (A19) includes boundary terms. Let $\mathsf{\Omega }$ be the spatial domain with boundary $\partial \mathsf{\Omega }$.

Appendix B.5.1. Diffusion Term Example

Consider the diffusion operator $L_d=-D{\nabla }^2$. Then:

$$\begin{array}{cc}\hfill D_E& =\langle \dot{E},-D{\nabla }^{2}E\rangle \hfill \\ \hfill & =-D{\int }_{\mathsf{\Omega }}\dot{E}{\nabla }^{2}EdV\hfill \\ \hfill & =D{\int }_{\mathsf{\Omega }}\nabla \dot{E}·\nabla EdV-D{\oint }_{\partial \mathsf{\Omega }}\dot{E}(\nabla E·\mathbf{n})dS\hfill \end{array}$$

(117)

Thus:

$$D_E=D_E^{\left(\mathrm{int}\right)}+{\mathsf{\Phi }}_{\mathrm{boundary}}$$

(A28)

where:

$$D_E^{\left(\mathrm{int}\right)}=D{\int }_{\mathsf{\Omega }}\nabla \dot{E}·\nabla EdV\ge 0$$

(A29)

is the internal dissipation, and:

$${\mathsf{\Phi }}_{\mathrm{boundary}}=-D{\oint }_{\partial \mathsf{\Omega }}\dot{E}(\nabla E·\mathbf{n})dS$$

(A30)

is the energy flux through the boundary.

Appendix B.5.2. General Statement

For any operator $L_d$ containing spatial derivatives:

$$\langle \dot{E},L_{d}E\rangle =D_E^{\left(\mathrm{int}\right)}+{\mathsf{\Phi }}_{\mathrm{boundary}}$$

(A31)

where $D_E^{\left(\mathrm{int}\right)}$ is the volumetric dissipation, and ${\mathsf{\Phi }}_{\mathrm{boundary}}$ is the energy flux through the boundary.

The energy continuity equation becomes:

$$\begin{array}{|c|}\hline {\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E^{\left(\mathrm{int}\right)}+{\mathsf{\Phi }}_{\mathrm{boundary}}=J_{S\to E}\\ \hline\end{array}$$

(A32)

Appendix B.6. Special Cases

Appendix B.6.1. Case 1: ODE Systems

For ODEs, there is no spatial domain, so ${\mathsf{\Phi }}_{\mathrm{boundary}}=0$. The equation reduces to:

$${\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E=J_{S\to E}$$

(A33)

where $D_E$ and $J_{S\to E}$ are simple products.

Appendix B.6.2. Case 2: No-Flux Boundaries

If the system has no-flux boundary conditions ($\nabla E·\mathbf{n}=0$ on $\partial \mathsf{\Omega }$), then ${\mathsf{\Phi }}_{\mathrm{boundary}}=0$, and the ODE form applies.

Appendix B.6.3. Case 3: Steady State

In steady state, $d{\mathcal{E}}_E/dt=0$, so:

$$D_E=J_{S\to E}$$

(A34)

All energy entering the environment from the source is dissipated.

Appendix B.7. The L 2 Norm Special Case

If we choose the specific energy functional:

$${\mathcal{E}}_E={\displaystyle \frac{1}{2}}{\parallel E\parallel }^2={\displaystyle \frac{1}{2}}{\int }_{\mathsf{\Omega }}{E}^{2}dV$$

(A35)

Then:

$${\displaystyle \frac{d{\mathcal{E}}_E}{dt}}=\langle E,\dot{E}\rangle$$

(A36)

and the energy continuity equation becomes:

$$\langle E,L_{c}E\rangle +\langle E,L_{d}E\rangle =\langle E,L_{d}S\rangle$$

(A37)

Appendix B.8. Connection to Thermodynamics

The energy continuity equation has a natural thermodynamic interpretation. For the environmental field:

• ${\mathcal{E}}_E$ is the accumulated environmental energy
• $D_E\ge 0$ is the entropy production rate
• $J_{S\to E}$ is the energy transfer rate from source to environment

The second law of thermodynamics requires $D_E\ge 0$ for passive systems. However, for active systems (e.g., the van der Pol oscillator), $D_E$ can be negative locally, representing energy injection.

Appendix B.9. Verification of the Derivation

We verify that Eq. (A24) is consistent with the fundamental equation $LE=L_{d}S$.

Proof. 

Starting from $LE=L_{d}S$:

$$\begin{array}{cc}\hfill \langle \dot{E},LE\rangle & =\langle \dot{E},L_{d}S\rangle \hfill \\ \hfill \langle \dot{E},L_{c}E\rangle +\langle \dot{E},L_{d}E\rangle & =J_{S\to E}\hfill \\ \hfill {\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E& =J_{S\to E}\hfill \end{array}$$

(A38)

which is exactly Eq. (A24). □

Appendix B.10. Summary

The energy continuity equation for the environmental field has been derived from first principles:

$$\begin{array}{|c|}\hline {\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E=J_{S\to E}\\ \hline\end{array}$$

(A39)

where:

$${\displaystyle \frac{d{\mathcal{E}}_E}{dt}}\equiv \langle \dot{E},L_{c}E\rangle ,D_E\equiv \langle \dot{E},L_{d}E\rangle ,J_{S\to E}\equiv \langle \dot{E},L_{d}S\rangle$$

(A40)

For PDEs, the dissipation term naturally splits into internal dissipation and boundary flux:

$$D_E=D_E^{\left(\mathrm{int}\right)}+{\mathsf{\Phi }}_{\mathrm{boundary}}$$

(A41)

yielding the complete form:

$$\begin{array}{|c|}\hline {\displaystyle \frac{d{\mathcal{E}}_E}{dt}}+D_E^{\left(\mathrm{int}\right)}+{\mathsf{\Phi }}_{\mathrm{boundary}}=J_{S\to E}\\ \hline\end{array}$$

(A42)

This completes the derivation of the energy continuity equation for the QSER framework.