R-matrix and uncertainties
==========================
The matrix :math:`R(s)` in observation space describes the error covariance of the model-observation comparison, assuming a fixed scaling factor :math:`s` in the model.
This section summarizes the assumptions made on uncertainties to construct :math:`R(s)`.
The construction of :math:`R(s)` is based on the work by :footcite:t:`Bruch2025`.

The uncertainty generally consists of multiple components.
We can approximate these as additive in the :math:`R(s)` matrix:

.. math::

   R^\text{combined}(s) = R^\text{transport}(s) + R^\text{flux}(s) + R^\text{bcens} + R^\text{diagonal} .

Here, :math:`R^\text{flux}(s)` captures the uncertainty distribution of fluxes within the flux categories used for the inversion (see :doc:`inversion`), :math:`R^\text{bcens}` is the ensemble-estimated uncertainty due to lateral boundary conditions, and :math:`R^\text{diagonal}` contains uncorrelated uncertainties such as observation uncertainties.

The construction of the R-matrix is implemented distributed in :py:class:`dubfi.fluxes.dataprovider_mpi_worker.MpiDistMecReaderWorker`, :py:class:`dubfi.linalg.mpi_worker.MpiDensePostRWorker`, :py:class:`dubfi.fluxes.dataprovider_dense.DataProvider`, :py:class:`dubfi.linalg.dense.DenseEnsLinPostR`, and :py:class:`dubfi.linalg.csc_sparse.CscEnsLinPostR`.
This construction depends on many tuning parameters, most of which are found in the configuration group :code:`uncertainty`.


Transport uncertainty
---------------------
We assume that we have a transport ensemble of :math:`M^\text{meteo}` members with known model predictions :math:`x^m_i(s)` that estimate the transport uncertainty.
In general, these model predictions will be approximations, which we will not discuss further here.
We further have a best model estimate :math:`x_i(s)`.

The first and main component for constructing :math:`R(s)` is the ensemble estimated transport uncertainty:

.. math::

   R^\text{transport}(s)_{ij} = \Lambda_{ij} \frac{1}{M^\text{meteo}-1} \sum_{m=1}^{M^\text{meteo}} \big[ x^m_i(s) - \bar{x}_i(s) \big] \big[ x^m_j(s) - \bar{x}_j(s) \big]

where :math:`\bar{x} = \frac{1}{M^\text{meteo}} \sum_m x^m` is the ensemble average and :math:`\Lambda` is a localization matrix.
We usually construct this localization as Gaussian suppression in horizontal, vertical, and temporal distance between observations (tuning parameters :code:`uncertainty.{hscale_m,vscale_m,tscale}`).
The transport uncertainty may be inflated by changing :code:`uncertainty.meteo_inflation`.


.. _approx_ensemble_members:

Approximating ensemble members
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Separating all flux categories in each ensemble member can be numerically very expensive.
We consider the cases where ensemble members group scaling factors to reduce the number of required tracers.
Consider a grouping matrix :math:`G_{kn}` where index :math:`k` selects a flux category (or scaling factor) and :math:`n` selects a group.
The columns of :math:`G` form basis vectors that are distinguished by the ensemble members.
We require that :math:`\sum_n G_{kn} = 1` for all :math:`k`.
We can then compute the concentration due to basis vector :math:`n` and denoted this by :math:`\hat{x}_n^m`.
To approximate :math:`x^m(s)`, we use the approximation (neglecting a possible additive background):

.. math::
   x^m(s) \approx \sum_{kn} G_{kn} s_k \frac{x^\text{det}_k}{\hat{x}_n^\text{det}} \hat{x}_n^m .

Here we use that the contribution :math:`x^\text{det}_k` of flux category :math:`k` in the deterministic run is known.
This approximation is used in the input data format required by DUBFI, see :ref:`input_format_ensemble`.
If you do not want to use this approximation, simply choose :math:`G` as the identity matrix.

Lateral boundary conditions
---------------------------
We consider an ensemble of :math:`M^\text{bc}` boundary conditions :math:`b^m_i` that represents the possible spread on the lateral boundaries. The resulting error covariance matrix is:

.. math::

   R^\text{bcens}_{ij} = \Lambda_{ij} \frac{1}{M^\text{bc}-1} \sum_{m=1}^{M^\text{bc}} \big[ b^m_i - \bar{b}_i \big] \big[ b^m_j - \bar{b}_j \big]

Flux ensemble
-------------
A flux ensemble is currently not implemented:

.. math::

   R^\text{flux} = 0

Uncorrelated uncertainties
--------------------------
Uncorrelated observation uncertainties only contribute to the diagonal of :math:`R(s)`.
We define

.. math::

   R^\text{diagonal}_{ij} = \delta_{ij} n_i^\text{nearby} \sigma_i^2, \qquad \sigma_i = \sigma_i^\text{obs} + \sigma_i^\text{baseline}

where :math:`\sigma_i^\text{baseline}` is a tuning parameter (:code:`uncertainty.uncorrelated`) that may depend on the station (:code:`uncertainty.uncorrelated_per_site`) and :math:`n_i^\text{nearby}` is the number of considered observations taken at the same time in very close distance.
The observation-based uncertainty :math:`\sigma_i^\text{obs}` is based on the temporal variability of the observed concentrations (:code:`uncertainty.obs_stddev_prefactor` and :code:`time_diff_prefactor`)

We set :math:`n_i^\text{nearby}` to the number of utilized sampling heights at the station and observation time of index :math:`i`.
For example, if we use simultaneous observations at three sampling heights of a tall tower observation site, then we weight the uncertainties with the factor :math:`\sqrt{3}` such that the combined uncertainty is comparable to the uncertainty of a single site.
This construction is based on the approximation that observations are precise and the uncertainty only needs to describe the limited ability of the model to reproduce the observations.

Uncertainty inflation
---------------------
Large deviations between model and observations can occur due to strong signals of wrong fluxes, or due to badly represented local effects.
To avoid a strong impact of such local effects on the inversion results, we need to ignore outliers or limit their influence by increasing their uncertainties.
Here we consider a fixed upper bound of the absolute model-observation mismatch (considering a priori, far-field-corrected model data) relative to the uncertainty.
Starting from :math:`R^\text{combined}(s)`, we introduce a weighting factor

.. math::

   \phi_i = \max\left\{1, \frac{|x_i(s^\text{prior}) - y_i|}{\kappa \sqrt{R^\text{combined}(s^\text{prior})_{ii}}} \right\}

where the tuning parameter :math:`\kappa` (:code:`data_filter.outlier_threshold`) is the maximum number of standard deviations by which the model may deviate from the observations.
In :math:`R^\text{weighted}(s)_{ij} = \phi_i \phi_j R^\text{combined}(s)`, each observation deviates from the prior by at most :math:`\kappa` standard deviations.

Inflation and final result
--------------------------
We inflate the observation uncertainties by a global factor :math:`\gamma` (:code:`uncertainty.global_inflation`) to reduce the overall weight of the uncertainties and obtain the final results :math:`R(s)_{ij} = \gamma^2 R^\text{weighted}(s)_{ij}`.

.. math::

   R(s)_{ij} = \gamma^2 \phi_i \phi_j \left[ R^\text{transport}(s) + R^\text{bcens} + R^\text{fluxens} + R^\text{diagonal} \right]_{ij}

Special case: Diagonal R matrix
-------------------------------
It is possible to discard the ensemble and use a diagonal :math:`R(s)` by choosing the implementation :code:`diagonal` in :py:mod:`dubfi.fluxes`.

.. footbibliography::