# An iterative optimization scheme to accommodate inequality constraints in air quality geostatistical estimation of multivariate PM

Maxime Beauchamp, Bertrand Bessagnet

► To cite this version:

Maxime Beauchamp, Bertrand Bessagnet. An iterative optimization scheme to accommodate inequality constraints in air quality geostatistical estimation of multivariate PM. Heliyon, 2023, pp.e17413. 10.1016/j.heliyon.2023.e17413. hal-04140956

HAL Id: hal-04140956

https://hal.science/hal-04140956

Submitted on 26 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Journal Pre-proof

An iterative optimization scheme to accommodate inequality constraints in air quality geostatistical estimation of multivariate PM

Maxime Beauchamp and Bertrand Bessagnet

Heliyon

PII: S2405-8440(23)04621-2

DOI: https://doi.org/10.1016/j.heliyon.2023.e17413

Reference: HLY e17413

To appear in: *Heliyon*

Received date: 22 August 2022

Revised date: 14 June 2023

Accepted date: 15 June 2023

Please cite this article as: M. Beauchamp and B. Bessagnet, An iterative optimization scheme to accommodate inequality constraints in air quality geostatistical estimation of multivariate PM, *Heliyon*, 9, e17413, doi: https://doi.org/10.1016/j.heliyon.2023.e17413.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2023 Published by Elsevier.
# An iterative optimization scheme to accommodate inequality constraints in air quality geostatistical estimation of multivariate PM

Maxime Beauchamp<sup>a,*,1</sup>, Bertrand Bessagnet<sup>b</sup>

## ARTICLE INFO

### Keywords:

- air quality

- cokriging

- optimization

- penalization

- multivariate particulate matter

- chemistry-transport model

## Abstract

The kriging-based estimation of the different types of atmospheric particulate matter (PM) pollutions defined in the air quality regulation raises some operational problems because the (co)kriging equations are obtained by minimizing a linear combination of the estimation variances subject to unbiasedness constraints. As a consequence, the estimation process can result in total PM₁₀ concentrations that are less than the PM₂.₅ concentrations which would be physically impossible. In a previous publication, it was shown that a convenient external drift modelling can reduce the number of spatial locations where the inequality constraint is not satisfied, without completely solving the problem. In this work, the formulation of the cokriging system is modified, inspired by previous works focusing on positive kriging. The introduction of additional constraints on the cokriging weights are presented, leading to a unique and optimal solution to the problem of cokriging under inequality constraints between two variables. Some computational and algorithmic details are introduced. An evaluation of the penalized cokriging is provided by using the European PM monitoring sites dataset: some maps and performance scores are given to assess the relevance of our iterative optimization scheme.

## 1. Introduction

In air quality, the geostatistical estimation is commonly used [9, 19, 20, 22, 30] to produce maps of ambient air regulatory pollutants [15], including PM₁₀ and PM₂.₅, the particles whose diameter are respectively smaller than 10 µg.m⁻³ and 2.5 µg.m⁻³. If the PM₁₀ mapping is a topic already well documented [4, 8], the interest for PM₂.₅ is growing [11, 34, 39] for public health reasons.

Measurements of PM₂.₅ having started later than PM₁₀ measurements, their monitoring network is less developed, which can affect the final mapping. In the French PREV'AIR system, see e.g. [33], the analyzed maps of ozone and PM₁₀ are built by kriging the observations [26], including simulations of the CTM CHIMERE model [25] as external drift. Applying such a technique to PM₂.₅ does not take into account the inequality between the concentrations of PM₂.₅ and PM₁₀, and can therefore lead to PM₂.₅ estimations greater than PM₁₀. Ad-hoc corrections were already considered, see e.g. [33] to estimate PM₂.₅ as the corresponding PM₁₀ times the ratio PM₂.₅/PM₁₀ simulated by CHIMERE. However, this method is not satisfactory given the uncertainties in the model-based simulations, in particular during air pollution episode. For the joint estimation of both pollutants, cokriging is the most common solution used with a linear coregionalization model [38]. Different solutions were already implemented: [35] proposed a Bayesian version of the linear model of coregionalization

applied to air quality data (CO, NO, NO₂). Multivariate non-stationarity processes were also addressed by Bayesian formulations, see e.g. [7], with spatial variations of the coefficients in the covariance model [18].

Beauchamp et al. [3] also developed a cokriging model for mapping PM₁₀ and PM₂.₅, in order to improve the precision of PM mapping. In this work, the local means of both variables are related by an additive model, but the latter does not ensure to satisfy the physical relationship existing between the two after their joint estimation. Indeed, the PM₂.₅ concentration is by definition less than PM₁₀. At most, the additive model of Beauchamp et al. [3] ensures the drift, a local regression of the observations on the Chemistry Transport Model (CTM) outputs [24], to be physically consistent thanks to the appropriate unbiasedness conditions. Thus, when the correlation between the model and the two observational datasets is good enough, the final estimation is very likely physically consistent. This is no longer the case if the drift poorly approaches the data. Anyway, this cokriging model with additive external drift has demonstrated its ability to better reproduce the PM concentration levels for mapping. It significantly reduces the number of cases in which PM₂.₅ estimations are greater than PM₁₀ (from 5 to 0.5%). However, even if the operational solution to deal with the inequality constraint is satisfactory, the mathematical problem associated is not entirely solved.

The way of dealing with inequalities in kriging has already been addressed, but mostly when the only information available is greater or less than a given value [16, 23]. For the estimation itself, Michalak [27] proposed a Gibbs sampler-

<figure><img src="image_2.png" /></figure>

<footnote><sup>a</sup>IMT Atlantique Bretagne-Pays de la Loire Campus de Brest Technopôle Brest-Iroise CS 83818 29238 Brest cedex 03 France</footnote>

<footnote><sup>b</sup>European Commission, Joint Research Centre (JRC), Ispra, Italy</footnote>

<footnote>*Maxime Beauchamp</footnote>

<footnote>✉ maxime.beauchamp@imt-atlantique.fr (M. Beauchamp); bertrand.bessagnet@ec.europa.eu (B. Bessagnet)</footnote>

<footnote>ORCID(s): 0000-0003-0605-8935 (M. Beauchamp)</footnote>
based approach for inequality-constrained geostatistical interpolation, supported by a review on the existing works about this topic: non-negativity of the weights [14], Langrange multipliers based approaches [28] or Monte-Carlo methodologies [1], etc. Additional equality and inequality constraints may also be used to solve kriging problems, such as negative weights in ordinary and simple kriging [37] or multiple indicator kriging [36].

In this work, a new algorithmic approach is presented to ensure PM<sub>2.5</sub> estimation be less than PM<sub>10</sub>, by introducing additional inequality constraints on the cokriging weights, writing the related optimization problem and solving it. The approach relies on Karush-Kuhn-Tucker conditions, inspired by the work of [2] on positive kriging. Additional constraints are formulated on the cokriging weights. When not met, the constraints result in rewriting the cokriging system. It yields an iterative procedure which is used until all the constraints are satisfied to ensure the consistency between both PM<sub>2.5</sub> and PM<sub>10</sub> estimation.

## 2. Methods

In Beauchamp et al. [3], the PM description is made through the simple additive model (1):

$$Z(\mathbf{x}) = Y(\mathbf{x}) + W(\mathbf{x}) \quad (1)$$

with $Y(\mathbf{x}) \le Z(\mathbf{x}) \forall \mathbf{x}$.

Here, the joint estimation between $Z(\mathbf{x})$, the PM<sub>10</sub> concentration and $Y(\mathbf{x})$, the PM<sub>2.5</sub> concentration, is investigated. Such a geostatistical estimation is known as cokriging [10].

The two components of the sum are supposed to be non-stationary, but can be explained by the use of deterministic covariates $f^i(\mathbf{x})$ and $g^j(\mathbf{x})$, see Eq. (2):

$$ 
\begin{aligned}
Y(\mathbf{x}) &= m_Y(\mathbf{x}) + R(\mathbf{x}) \\
&= a_0 + \sum_{j=1}^{l} a_j f^i(\mathbf{x}) + R(\mathbf{x}) \\
W(\mathbf{x}) &= m_W(\mathbf{x}) + T(\mathbf{x}) \\
&= b_0 + \sum_{j=1}^{p} a_j g^j(\mathbf{x}) + T(\mathbf{x})
\end{aligned}
\quad (2) $$

$f^i$, $i = 1, \dots, l$ and $g^j$, $j = 1, \dots, p$ are the covariates respectively used in the computation of the local means $m_Y(\mathbf{x})$ and $m_W(\mathbf{x})$. $l$ and $p$ are the number of covariates for the regression of $Y$ on the $f^i$ and of $W$ on the $g^j$. Last, $R$ and $T$ are residuals assumed to be second-order stationary random functions with zero mean. This framework can be used for any similar modelling, and even simplified if not using any covariates.

Let denote $Y^{\text{CK}}(\mathbf{x}_0)$ and $Z^{\text{CK}}(\mathbf{x}_0)$ the cokriging-based estimation of PM<sub>2.5</sub> and PM<sub>10</sub> at location $\mathbf{x}_0$. For sake of simplicity, $Y$ and $Z$ respectively refer to variables with index 1 and 2 throughout the paper. We aim at producing a joint estimation such that $Y^{\text{CK}}(\mathbf{x}_0)$ is less than $Z^{\text{CK}}(\mathbf{x}_0)$. The additive relationship between $Z$ and $Y$ implies that direct and cross covariances are linked through the following set of equations:

$$ 
\begin{cases}
C_Z(\mathbf{h}) = C_R(\mathbf{h}) \\
C_{Z,Y}(\mathbf{h}) = C_R(\mathbf{h}) + C_{R,T}(\mathbf{h}) \\
C_{Y,Z}(\mathbf{h}) = C_R(\mathbf{h}) + C_{T,R}(\mathbf{h}) \\
C_Y(\mathbf{h}) = C_R(\mathbf{h}) + C_{R,T}(\mathbf{h}) \\
\qquad + C_{R,T}(-\mathbf{h}) + C_T(\mathbf{h})
\end{cases}
$$

The same conditions can also be written using variograms:

$$ 
\begin{cases}
\gamma_Z(\mathbf{x}, \mathbf{x} + \mathbf{h}) = C_R(0) - C_R(\mathbf{h}) \\
\gamma_{Y,Z}(\mathbf{x}, \mathbf{x} + \mathbf{h}) = \gamma_R(\mathbf{h}) + \gamma_{R,T}(\mathbf{h}) \\
\gamma_Y(\mathbf{x}, \mathbf{x} + \mathbf{h}) = \gamma_R(\mathbf{h}) + \gamma_T(\mathbf{h}) \\
\qquad + 2\gamma_{R,T}(\mathbf{h})
\end{cases}
$$

$C(\cdot)$ denotes a covariance and $\gamma(\cdot)$ a variogram, see e.g. Cressie and Wikl [13]. In Beauchamp et al. [3], it is described why the use of variograms instead of covariances is preferable and well suited to the data and their related hypothesis.

$Z^{\text{CK}}(\mathbf{x}_0)$ is obtained by solving the cokriging system (CK), see Eq. (3) and Sect. 2.5 in [3]:

$$Z^{\text{CK}}(\mathbf{x}_0) = \sum_{\beta=1}^{M} \nu_{\beta}' Z(\mathbf{x}_{\beta}) + \sum_{\alpha=1}^{N} \lambda_{\alpha}' Y(\mathbf{x}_{\alpha}) \quad (3)$$

where $\alpha = 1, \dots, N$ and $\beta = 1, \dots, M$ are the indices for the observations $Y(\mathbf{x}_{\alpha})$ and $Z(\mathbf{x}_{\beta})$. $\nu_{\beta}' = \nu'(\mathbf{x}_{\beta})$ and $\lambda_{\alpha}' = \lambda'(\mathbf{x}_{\alpha})$ are the cokriging weights of the observations $Z(\mathbf{x}_{\beta})$ and $Y(\mathbf{x}_{\alpha})$. Modifying a bit the cokriging system notations of Beauchamp et al. [3] to ease the link with this work, see Eq. (4), the cokriging system of $Z^{\text{CK}}(\mathbf{x}_0)$ becomes:

$$ 
\begin{bmatrix}
\mathbf{K}_{22} & \mathbf{K}_{21} & 1 & \mathbf{g}_2^j & 0 & \mathbf{f}_2^i \\
\mathbf{K}_{12} & \mathbf{K}_{11} & 0 & 0 & 1 & \mathbf{f}_1^i \\
1 & 0 & 0 & 0 & 0 & 0 \\
\mathbf{g}_2^j & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 \\
\mathbf{f}_2^i & \mathbf{f}_1^i & 0 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\nu' \\
\lambda' \\
\mu_0 \\
\mu_j \\
\nu_0 \\
\nu_i
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{K}_{20} \\
\mathbf{K}_{120} \\
1 \\
g_0^j \\
0 \\
f_0^i
\end{bmatrix}
\quad (4) $$

where:

$$ \mathbf{K}_{11} = \{C_Y(\mathbf{x}_{\alpha} - \mathbf{x}_{\alpha'})\}. $$

$$ \mathbf{K}_{22} = \{C_Z(\mathbf{x}_{\beta} - \mathbf{x}_{\beta'})\}. $$

$$ \mathbf{K}_{12} = \{C_{Y,Z}(\mathbf{x}_{\alpha} - \mathbf{x}_{\beta})\}. $$

$$ \mathbf{K}_{20} = \{C_Z(\mathbf{x}_{\beta} - \mathbf{x}_{0})\}. $$
$$ \mathbf{K}_{120} = \{C_{Y,Z}(\mathbf{x}_{\alpha} - \mathbf{x}_{0})\}. $$

$v' = (v'_{\beta})$ and $\lambda' = (\lambda'_{\alpha})$ are the cokriging weights of the observations $Z(\mathbf{x}_{\beta})$ and $Y(\mathbf{x}_{\alpha})$. $\mathbf{g}_{2}^{j}$ denotes the vector of the j$^{th}$ covariate $g^{j}$ at locations $\mathbf{x}_{\beta}$ {$g^{j}(\mathbf{x}_{\beta})$}. $\mathbf{f}_{2}^{i}$ denotes the vector of the i$^{th}$ covariate $f^{i}$ at locations $\mathbf{x}_{\beta}$ {$f^{i}(\mathbf{x}_{\beta})$}. $\mathbf{f}_{1}^{i}$ denotes the vector of the i$^{th}$ covariate $f^{i}$ at locations $\mathbf{x}_{\alpha}$ {$f^{i}(\mathbf{x}_{\alpha})$}. $g_{0}^{j}$ and $f_{0}^{i}$ respectively denote the j$^{th}$ and i$^{th}$ covariates $g^{j}$ and $f^{i}$ at the target location $\mathbf{x}_{0}$. Last, $\mu_{0}, \mu_{j}, v_{0}, v_{i}$ are the Lagrange multipliers for the unbiasedness conditions, see Eq.(11a) to Eq.(11d) in Beauchamp et al. [3].

For sake of generality, this system is given with covariances but can be easily written with variograms, as this tool was used in Beauchamp et al. [3] and again in this study. As denoted by Isaaks and Srivastava [21], this type of cokriging estimation may lead to negative values and the effect of the secondary variable on the estimation is weak depending on the spatial sampling of the two variables, but it is not relevant here because there is generally more PM$_{10}$ data available than PM$_{2.5}$ and the deterministic covariate almost ensures the positivity of the estimation, see again Beauchamp et al. [3].

In a more convenient synthetic notations:

$$ \begin{aligned} & \mathbf{KP} + \mathbf{AM} = \mathbf{K}_{0} \\ & \mathbf{A}'\mathbf{P} = \mathbf{F}_{0} \end{aligned} \quad (5) $$

where:

$$ \mathbf{K} = \begin{bmatrix} \mathbf{K}_{11} & \mathbf{K}_{12} \\ \mathbf{K}_{21} & \mathbf{K}_{22} \end{bmatrix}, \quad \mathbf{A} = \begin{bmatrix} 1 & \mathbf{g}_{2}^{j} & 0 & \mathbf{f}_{2}^{i} \\ 0 & 0 & 1 & \mathbf{f}_{1}^{i} \end{bmatrix}, $$

$$ \mathbf{P} = \begin{bmatrix} v'_{\beta} \\ \lambda'_{\alpha} \end{bmatrix}, \quad \mathbf{M} = \begin{bmatrix} \mu_{0} \\ \mu_{j} \\ v_{0} \\ v_{i} \end{bmatrix} $$

$$ \mathbf{K}_{0} = \begin{bmatrix} \mathbf{K}_{20} \\ \mathbf{K}_{120} \end{bmatrix}, \quad \mathbf{F}_{0} = \begin{bmatrix} 1 \\ g_{0}^{j} \\ 0 \\ f_{0}^{i} \end{bmatrix} $$

Let suppose now that $Z^{\text{CK}}(\mathbf{x}_{0})$ from Eq. (5) is already computed. The estimation $Y^{\text{CK}}(\mathbf{x}_{0})$ writes:

$$ Y^{\text{CK}}(\mathbf{x}_{0}) = \sum_{\alpha=1}^{N} \lambda_{\alpha} Y(\mathbf{x}_{\alpha}) + \sum_{\beta=1}^{M} v_{\beta} Z(\mathbf{x}_{\beta}) $$

where the weights $\lambda_{\alpha}$ of $Y(\mathbf{x}_{\alpha})$ and the weights $v_{\beta}$ of $Z(\mathbf{x}_{\beta})$ are found by solving a cokriging system similar to (4), wherein the order of the variables is simply modified accordingly in

matrix $\mathbf{K}$ and vectors $\mathbf{P}, \mathbf{M}$ and $\mathbf{K}_{0}$. The unbiasedness conditions also modify the contents of matrix $\mathbf{A}$:

$$ \mathbf{A} = \begin{bmatrix} 1 & 0 & 0 & \mathbf{f}_{1}^{i} \\ 0 & \mathbf{g}_{2}^{j} & 1 & \mathbf{f}_{2}^{i} \\ 0 & 0 & 1 & f_{0}^{i} \end{bmatrix}, \quad \mathbf{F}_{0} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ f_{0}^{i} \end{bmatrix} $$

## 2.1. Additional constraints on the weights

$Y^{\text{CK}}(\mathbf{x}_{0})$ has to be less than the prior estimation $Z^{\text{CK}}(\mathbf{x}_{0})$ of the first variable, i.e.:

$$ \begin{aligned} Y^{\text{CK}}(\mathbf{x}_{0}) &\le Z^{\text{CK}}(\mathbf{x}_{0}) \\ &\le \sum_{\beta=1}^{M} v'_{\beta} Z(\mathbf{x}_{\beta}) + \sum_{\alpha=1}^{N} \lambda'_{\alpha} Y(\mathbf{x}_{\alpha}) \end{aligned} \quad (6) $$

Let precise that cutting the estimation of $Y^{\text{CK}}$ as the minimum between $(Y^{\text{CK}}, Z^{\text{CK}})$ is not a good solution because even if it can lead to satisfying results, both in terms of mapping and obviously being consistent with the physical constraints, it does not enable to compute the corresponding cokriging standard deviation.

Replacing the cokriging of $Z$ by its mathematical expression given in Eq.3 leads to:

$$ \sum_{\alpha=1}^{N} \lambda_{\alpha} Y(\mathbf{x}_{\alpha}) \le \sum_{\alpha=1}^{N} \lambda'_{\alpha} Y(\mathbf{x}_{\alpha}) + \sum_{\beta=1}^{M} (v'_{\beta} - v_{\beta}) Z(\mathbf{x}_{\beta}) $$

Revisiting this inequality, we make the choice of evenly decomposing $\sum_{\beta=1}^{M} (v'_{\beta} - v_{\beta}) Z(\mathbf{x}_{\beta}) = \sum_{\alpha=1}^{N} \tilde{\lambda}'_{\alpha}$ as the sum of $N$ terms where:

$$ \tilde{\lambda}'_{\alpha} = \left[ \sum_{\beta=1}^{M} (v'_{\beta} - v_{\beta}) Z(\mathbf{x}_{\beta}) \right] \frac{l(\mathbf{x}_{\alpha})}{Y(\mathbf{x}_{\alpha})}, \quad \alpha = 1, \dots, N $$

with $\sum_{\alpha=1}^{N} l(\mathbf{x}_{\alpha}) = 1$ and $Y(\mathbf{x}_{\alpha}) > 0 \forall \alpha$. In this work, we consider a uniform repartition scheme, then $l(\mathbf{x}_{\alpha}) = 1/N \forall \mathbf{x}_{\alpha}$, but any other scheme could be experimented. The additional weights $\tilde{\lambda}'_{\alpha}$ are then distributed over the $N$ weights $\lambda'_{\alpha}$:

$$ \sum_{\alpha=1}^{N} \lambda_{\alpha} Y(\mathbf{x}_{\alpha}) \le \sum_{\alpha=1}^{N} [\lambda'_{\alpha} + \tilde{\lambda}'_{\alpha}] Y(\mathbf{x}_{\alpha}) $$

Thus, at $v_{\beta}$ fixed, the condition (7):

$$ \lambda_{\alpha} \le \lambda'_{\alpha} + \tilde{\lambda}'_{\alpha} \quad \forall \alpha \quad (7) $$

is sufficient to ensure the estimation $Y^{\text{CK}}(\mathbf{x}_{0})$ to be less than $Z^{\text{CK}}(\mathbf{x}_{0})$.
Because of both unbiasedness cokriging conditions

$$ \sum_{\alpha=1}^{N} \lambda_{\alpha} = 1 \text{ and } \sum_{\alpha=1}^{N} \lambda'_{\alpha} = 0, \text{ we have:} \\ \sum_{\alpha=1}^{N} \tilde{\lambda}'_{\alpha} \geq 1 + \lambda_{0} \quad \text{with } \lambda_{0} \geq 0. \qquad (8) $$

From the expression of $\tilde{\lambda}'_\alpha$ and inequality (8), we can use a similar framework to split and distribute a given quantity over the weights $v'_\beta$, we have:

$$ \sum_{\beta=1}^{M} v_{\beta} Z(\mathbf{x}_{\beta}) \leq \sum_{\beta=1}^{M} [v'_{\beta} - \tilde{v}'_{\beta}] Z(\mathbf{x}_{\beta}) \qquad (9) $$

with:

$$ \tilde{v}'_{\beta} = \frac{k(\mathbf{x}_{\beta})(1 + \lambda_{0})}{Z(\mathbf{x}_{\beta}) \sum_{\alpha=1}^{N} l(\mathbf{x}_{\alpha})/Y(\mathbf{x}_{\alpha})}, $$

and $\sum_{\beta=1}^{M} k(\mathbf{x}_\beta) = 1$. Again, if a uniform scheme distribution is used over the weights $v'_\beta$, then $k(\mathbf{x}_\beta) = 1/M \forall \mathbf{x}_\beta$.

A sufficient condition to satisfy inequality (9) is:

$$ v_{\beta} \leq v'_{\beta} - \tilde{v}'_{\beta} \quad \forall \mathbf{x}_{\beta} \qquad (10) $$

And because of condition (10) and both unbiasedness conditions

$$ \sum_{\beta=1}^{M} v_{\beta} = 0 \text{ and } \sum_{\beta=1}^{M} v'_{\beta} = 1: \\ \sum_{\beta=1}^{M} \tilde{v}'_{\beta} \leq 1 - v_{0} \quad \text{with } v_{0} \geq 0 \qquad (11) $$

Thus, for a given pair $(\lambda_0, v_0)$ satisfying conditions (8) and (11), the set of weights $\lambda_\alpha$ and $v_\beta$ compliant with inequality (6) have the following property:

$$ \frac{1 + \lambda_0}{\sum_{\alpha} l(\mathbf{x}_{\alpha})/Y(\mathbf{x}_{\alpha})} \leq \frac{1 - v_0}{\sum_{\beta} k(\mathbf{x}_{\beta})/Z(\mathbf{x}_{\beta})} \qquad (12) $$

Because the variables $Y(\mathbf{x})$ and $Z(\mathbf{x})$ are always positive, it is in most cases possible to find a subset of the observational data $\{Y(\mathbf{x}_\alpha), Z(\mathbf{x}_\beta)\}$ satisfying condition (12).

Regarding the other unbiasedness condition, Eqs. (13a,13b,13c,13d)

appearing in the two cokriging estimations, i.e.:

$$ \left\{ \begin{array}{l} \sum_{\alpha} \lambda_{\alpha} f^{i}(\mathbf{x}_{\alpha}) + \sum_{\beta} v_{\beta} f^{i}(\mathbf{x}_{\beta}) = f^{i}(\mathbf{x}_{0}) \quad (13a) \\ \\ \sum_{\alpha} \lambda'_{\alpha} f^{i}(\mathbf{x}_{\alpha}) + \sum_{\beta} v'_{\beta} f^{i}(\mathbf{x}_{\beta}) = f^{i}(\mathbf{x}_{0}) \quad (13b) \\ \\ \sum_{\beta} v_{\beta} g^{j}(\mathbf{x}_{\beta}) = 0 \quad (13c) \\ \\ \sum_{\beta} v'_{\beta} g^{j}(\mathbf{x}_{\beta}) = g^{j}(\mathbf{x}_{0}) \quad (13d) \end{array} \right. $$

In what follows, the notations $\overline{\lambda}_\alpha$, Eq. (14a), and $\overline{v}_\beta$, Eq. (14b), denote the penalized conditions (7) and (10):

$$ \overline{\lambda}_{\alpha} = \lambda'_{\alpha} + \tilde{\lambda}'_{\alpha} \qquad (14a) $$

$$ \overline{v}_{\beta} = v'_{\beta} - \tilde{v}'_{\beta} \qquad (14b) $$

## 2.2. The optimization problem

Here, we present the optimization problem (P) derived from the formulation of cokriging with inequality constraints:

$$ \left\{ \begin{array}{l} \Omega = (\lambda, v) \\ = \operatorname{argmin} \operatorname{Var}[Y^{\mathrm{CK}}(\mathbf{x}_0) - Y(\mathbf{x}_0)] | h_i(\Omega) = 0, g_j(\Omega) \le 0 \\ \text{where } h_i(\Omega) = \begin{cases} \sum_{\alpha} \lambda_{\alpha} - 1 = 0 \\ \sum_{\beta} v_{\beta} = 0 \\ \sum_{\alpha} \lambda_{\alpha} f^{i}(\mathbf{x}_{\alpha}) + \sum_{\beta} v_{\beta} f^{i}(\mathbf{x}_{\beta}) = f_{0}^{i} \\ \sum_{\beta} v_{\beta} g^{j}(\mathbf{x}_{\beta}) = 0 \end{cases} \\ \text{and } g_j(\Omega) = \begin{cases} \lambda_{\alpha} - \overline{\lambda}_{\alpha} \le 0, \quad \alpha = 1, \dots, N \\ v_{\beta} - \overline{v}_{\beta} \le 0, \quad \beta = 1, \dots, M \end{cases} \end{array} \right. $$

The Lagrangian of problem (P) writes:

$$ \begin{aligned} L(\Omega, \mu, \tau) &= \operatorname{Var}[Y^{\mathrm{CK}}(\mathbf{x}_0) - Y(\mathbf{x}_0)] + \sum_{i=1}^{l+p+2} \kappa_i h_i(\Omega) \\ &\quad + \sum_{j=1}^{2} \eta_j g_j(\Omega) \\ &= \operatorname{Var}[Y^{\mathrm{CK}}(\mathbf{x}_0) - Y(\mathbf{x}_0)] + 2\mu_0 \left( \sum_{\alpha} \lambda_{\alpha} - 1 \right) \\ &\quad + 2v_0 \sum_{\beta} (v_{\beta}) + 2\mu_j \left( \sum_{\beta} v_{\beta} g^{j}(\mathbf{x}_{\beta}) \right) \\ &\quad + 2v_i \left( \sum_{\alpha} \lambda_{\alpha} f^{i}(\mathbf{x}_{\alpha}) + \sum_{\beta} v_{\beta} f^{i}(\mathbf{x}_{\beta}) - f_{0}^{i} \right) \\ &\quad + \sum_{\alpha} \tau_{\alpha} (\lambda_{\alpha} - \overline{\lambda}_{\alpha}) + \sum_{\beta} \omega_{\beta} (v_{\beta} - \overline{v}_{\beta}) \end{aligned} $$
If there is a local minimum $\Omega^*$ of $\text{Var}[Y(\mathbf{x}_0) - Y^{\text{CK}}(\mathbf{x}_0)]$, it exists $\kappa = (\mu, \nu)$, $\mu \in \mathbb{R}^{l+1}$, $\nu \in \mathbb{R}^{p+1}$, and $\eta = (\tau, \omega)$, $\tau \in \mathbb{R}^N$, and $\omega \in \mathbb{R}^M$ so that:

$$
\left\{
\begin{aligned}
\nabla_{\Omega} L(\Omega^*, \mu, \tau) &= \begin{bmatrix} \frac{\partial L}{\partial \lambda^\alpha} \\ \frac{\partial L}{\partial \nu^\beta} \end{bmatrix} \\
&= \begin{bmatrix} \sum_{\alpha'} \lambda_{\alpha'} C_{\alpha\alpha'}^1 - C_{\alpha\alpha}^1 + \sum_{\beta} \nu_{\beta} C_{\alpha\beta}^{12} + \mu_0 + \nu_i f^i(\mathbf{x}_\alpha) + \tau_\alpha \\ \sum_{\beta'} \nu_{\beta'} C_{\beta\beta'}^2 - C_{\beta x}^{12} + \sum_\alpha \lambda_\alpha C_{\alpha\beta}^{12} + \nu_0 + \mu_j g^j(\mathbf{x}_\beta) + \nu_i f^i(\mathbf{x}_\beta) + \omega_\beta \end{bmatrix} \\
&= 0 \\
h_i(\Omega^*) = 0, \text{ i.e.} \quad & \begin{cases} \sum_\alpha \lambda_\alpha - 1 = 0 \\ \sum_\beta \nu_\beta = 0 \\ \sum_\alpha \lambda_\alpha f^i(\mathbf{x}_\alpha) + \sum_\beta \nu_\beta f^i(\mathbf{x}_\beta) = f_0^i \\ \sum_\beta \nu_\beta g^j(\mathbf{x}_\beta) = 0 \end{cases} \\
\eta_j g_j(\Omega^*) = 0, \text{ i.e.} \quad & \begin{cases} \forall \alpha, \tau_\alpha(\lambda_\alpha - \bar{\lambda}_\alpha) = 0, \text{ i.e. } \sum_\alpha \tau_\alpha(\lambda_\alpha - \bar{\lambda}_\alpha) = 0 \\ \forall \beta, \omega_\beta(\nu_\beta - \bar{\nu}_\beta) = 0, \text{ i.e. } \sum_\beta \omega_\beta(\nu_\beta - \bar{\nu}_\beta) = 0 \end{cases} \\
\tau_\alpha \ge 0, \alpha = 1, \dots, N \text{ and } \omega_\beta \ge 0, \beta = 1, \dots, M
\end{aligned}
\right.
$$

For sake of simplicity and generality, index 1 and 2 respectively denotes variable $Y$ and $Z$, and the notation $C_{\alpha\alpha'}^1$ stands for $C_1(|\mathbf{x}_\alpha - \mathbf{x}_{\alpha'}|)$ in what follows. The same applies for the other covariances.

These are the so-called Karush-Kuhn-Tucker (KKT) conditions, see e.g. Rothenberg [32], that can be written in a more readable matrix form as follows:

$$ \begin{bmatrix} \mathbf{K} & \mathbf{A} & \mathbf{I} & \mathbf{0} \\ \mathbf{A}' & 0 & 0 & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{P} \\ \mathbf{M} \\ \tau_{\alpha} \\ \omega_{\beta} \end{bmatrix} = \begin{bmatrix} \mathbf{K}_{0} \\ \mathbf{F}_{0} \end{bmatrix} \quad (15a) $$

where notations $\mathbf{K}$, $\mathbf{A}$, $\mathbf{P}$, $\mathbf{M}$, $\mathbf{K}_0$ and $\mathbf{F}_0$ are taken from Eq.(5) and with side constraints :

$$ \begin{align} \sum_{\alpha} \tau_{\alpha}(\lambda_{\alpha} - \bar{\lambda}_{\alpha}) &= 0 \tag{15b} \\ \sum_{\beta} \omega_{\beta}(\nu_{\beta} - \bar{\nu}_{\beta}) &= 0 \notag \end{align} $$

$$ \lambda_{\alpha} - \bar{\lambda}_{\alpha} \le 0 \quad \text{and} \quad \tau_{\alpha} \ge 0, \alpha = 1, \dots, N \\ \nu_{\beta} - \bar{\nu}_{\beta} \le 0 \quad \omega_{\beta} \ge 0, \beta = 1, \dots, M \quad (15c) $$

This estimation of $Y(\mathbf{x}_0)$ obtained by solving this KKT-based penalization, Eqs. (15b) and (15c), of the cokriging system (15a) is now denoted $Y^{\text{PCK}}(\mathbf{x}_0)$.

## 2.3. The complementary slackness conditions

Since the conditions (15c) requires all of the $\tau_\alpha$ and $\omega_\beta$ be zero or positive while all the terms $(\lambda_\alpha - \bar{\lambda}_\alpha)$ and $(\nu_\beta - \bar{\nu}_\beta)$ be zero or negative, the condition (15b) requires :

$$ \tau_{\alpha} = 0 \text{ or } \lambda_{\alpha} = \bar{\lambda}_{\alpha} \\ \text{and } \omega_{\beta} = 0 \text{ or } \nu_{\beta} = \bar{\nu}_{\beta} $$

This is a particular case of complementary slackness condition, see e.g. Boyd and Vandenberghe [6], where the variables $\tau_\alpha$ and $\omega_\beta$ are respectively the complementary pairs, also called Lagrangian multipliers or dual variables, corresponding to the constraints $\lambda_\alpha \le \bar{\lambda}_\alpha$ and $\nu_\beta \le \bar{\nu}_\beta$. Specifically, $\sum_\alpha \tau_\alpha(\lambda_\alpha - \bar{\lambda}_\alpha) = 0$ indicates that, at the optimal solution $\lambda_\alpha$ and without loss of generality, either $\tau_\alpha$ is zero or $\lambda_\alpha = \bar{\lambda}_\alpha$, i.e. the inequality constraint is binding. The second constraint $\sum_\beta \omega_\beta(\nu_\beta - \bar{\nu}_\beta) = 0$ can be interpreted in a similar way.

In what follows, the positive kriging approach of Barnes and Johnson [2] is adapted to address our problem: for a given solution (not necessarily optimal), the cokriging weights satisfying condition (15c) are denoted as the **basic weights $\lambda_\alpha$ and $\nu_\beta$**:

$$ \begin{align} B_{\alpha} &= \{\alpha, \lambda_{\alpha} < \bar{\lambda}_{\alpha}, \tau_{\alpha} = 0\} \notag \\ B_{\beta} &= \{\beta, \nu_{\beta} < \bar{\nu}_{\beta}, \omega_{\beta} = 0\} \tag{16} \end{align} $$
while the weights equal to their authorized upper-boundaries are identified as the **non-basic** weights:

$$ 
\begin{aligned}
\bar{B}_{\alpha} &= \{\alpha, \lambda_{\alpha} = \bar{\lambda}_{\alpha}, \tau_{\alpha} \ge 0\} \\
\bar{B}_{\beta} &= \{\beta, \nu_{\beta} = \bar{\nu}_{\beta}, \omega_{\beta} \ge 0\}
\end{aligned}
\qquad (17) $$

Along this line, the set of $\tau_\alpha$ and $\omega_\beta$ respectively complementary with the "basic" $\lambda_\alpha$ and $\nu_\beta$ will be called the "basic" $\tau_\alpha$ and $\omega_\beta$. The same idea applies for the definition of the non-basic $\tau_\alpha$ and $\omega_\beta$.

The penalized cokriging system for the estimation $Y^{\text{PCK}}(\mathbf{x}_0)$ can be rearranged and partitioned by ordering rows and columns according to the basic and non-basic components of the solution:

$$ \lambda = \begin{bmatrix} \lambda_b \\ \lambda_n \end{bmatrix}, \quad \nu = \begin{bmatrix} \nu_b \\ \nu_n \end{bmatrix}, \quad \tau = \begin{bmatrix} \tau_b \\ \tau_n \end{bmatrix}, \quad \omega = \begin{bmatrix} \omega_b \\ \omega_n \end{bmatrix} $$

where $\lambda_b$ denotes the set of basic weights $\lambda_\alpha$, $\alpha \in \mathcal{B}_\alpha$ and $\lambda_n$, $\alpha \in \bar{\mathcal{B}}_\alpha$ the set of non-basic weights $\lambda_\alpha$. The same type of notations are used for $\nu_\beta$, $\tau_\alpha$ and $\omega_\beta$.

Reordering rows and columns with the basic and non-basic nomenclature, Eqs. (16) and (17), gives:

$$ 
\begin{aligned}
& \begin{bmatrix} \mathbf{K}_{11}^{\text{bb}} & \mathbf{K}_{11}^{\text{bn}} & \mathbf{K}_{12}^{\text{bb}} & \mathbf{K}_{12}^{\text{bn}} \\ \mathbf{K}_{11}^{\text{nb}} & \mathbf{K}_{11}^{\text{nn}} & \mathbf{K}_{12}^{\text{nb}} & \mathbf{K}_{12}^{\text{nn}} \\ \mathbf{K}_{21}^{\text{bb}} & \mathbf{K}_{21}^{\text{bn}} & \mathbf{K}_{22}^{\text{bb}} & \mathbf{K}_{22}^{\text{bn}} \\ \mathbf{K}_{21}^{\text{nb}} & \mathbf{K}_{21}^{\text{nn}} & \mathbf{K}_{22}^{\text{nb}} & \mathbf{K}_{22}^{\text{nn}} \end{bmatrix} \begin{bmatrix} \lambda_b \\ \lambda_n \\ \nu_b \\ \nu_n \end{bmatrix} \\
& + \begin{bmatrix} 1 & 0 & 0 & \mathbf{f}_{1b}^i \\ 1 & 0 & 0 & \mathbf{f}_{1n}^i \\ 0 & \mathbf{g}_{2b}^j & 1 & \mathbf{f}_{2b}^i \\ 0 & \mathbf{g}_{2n}^j & 1 & \mathbf{f}_{2n}^i \end{bmatrix} \begin{bmatrix} \mu_0 \\ \mu_j \\ \nu_0 \\ \nu_i \end{bmatrix} \\
& + \begin{bmatrix} \mathbf{1}_{11}^{\text{bb}} & \mathbf{0}_{11}^{\text{bn}} & \mathbf{0}_{12}^{\text{bb}} & \mathbf{0}_{12}^{\text{bn}} \\ \mathbf{0}_{11}^{\text{nb}} & \mathbf{1}_{11}^{\text{nn}} & \mathbf{0}_{12}^{\text{nb}} & \mathbf{0}_{12}^{\text{nn}} \\ \mathbf{0}_{21}^{\text{bb}} & \mathbf{0}_{21}^{\text{bn}} & \mathbf{1}_{22}^{\text{bb}} & \mathbf{0}_{22}^{\text{bn}} \\ \mathbf{0}_{21}^{\text{nb}} & \mathbf{0}_{21}^{\text{nn}} & \mathbf{0}_{22}^{\text{nb}} & \mathbf{1}_{22}^{\text{nn}} \end{bmatrix} \begin{bmatrix} \tau_b \\ \tau_n \\ \omega_b \\ \omega_n \end{bmatrix} = \begin{bmatrix} \mathbf{K}_{10}^b \\ \mathbf{K}_{10}^n \\ \mathbf{K}_{120}^b \\ \mathbf{K}_{120}^n \end{bmatrix} \quad (18a)
\end{aligned}
$$

$$ \begin{bmatrix} 1 & 1 & 0 & 0 \\ 0 & 0 & \mathbf{g}_{2b}^j & \mathbf{g}_{2n}^j \\ 0 & 0 & 1 & 1 \\ \mathbf{f}_{1b}^i & \mathbf{f}_{1n}^i & \mathbf{f}_{2b}^i & \mathbf{f}_{2n}^i \end{bmatrix} \begin{bmatrix} \lambda_b \\ \lambda_n \\ \nu_b \\ \nu_n \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ f_0^i \end{bmatrix} \quad (18b) $$

$$ \begin{bmatrix} \tau_b & \tau_n & \omega_b & \omega_n \end{bmatrix} \begin{bmatrix} \lambda_b - \bar{\lambda}_b \\ \lambda_n - \bar{\lambda}_n \\ \nu_b - \bar{\nu}_b \\ \nu_n - \bar{\nu}_n \end{bmatrix} = 0 \quad (18c) $$

$$ \lambda_{\alpha} - \bar{\lambda}_{\alpha} \le 0 \quad \tau_{\alpha}^* \ge 0, \quad \alpha = 1, \dots, N \\ \nu_{\beta} - \bar{\nu}_{\beta} \le 0 \quad \omega_{\beta}^* \ge 0, \quad \beta = 1, \dots, M \quad (18d) $$

$$ \mathbf{K}_{ij}^{\text{bb}} = \{C_{ij}(\mathbf{x}_{\alpha} - \mathbf{x}_{\alpha}')\}, \alpha \in \mathcal{B}_{\alpha}, \alpha' \in \bar{\mathcal{B}}_{\alpha}. $$

$$ \mathbf{K}_{ij}^{\text{bn}} = \{C_{ij}(\mathbf{x}_{\alpha} - \mathbf{x}_{\alpha}')\}, \alpha \in \bar{\mathcal{B}}_{\alpha}, \alpha' \in \mathcal{B}_{\alpha}. $$

$$ \mathbf{K}_{ij}^{\text{nb}} = \{C_{ij}(\mathbf{x}_{\alpha} - \mathbf{x}_{\alpha}')\}, \alpha \in \mathcal{B}_{\alpha}, \alpha' \in \bar{\mathcal{B}}_{\alpha}. $$

$$ \mathbf{K}_{ij}^{\text{nn}} = \{C_{ij}(\mathbf{x}_{\alpha} - \mathbf{x}_{\alpha}')\}, \alpha \in \bar{\mathcal{B}}_{\alpha}, \alpha' \in \mathcal{B}_{\alpha}. $$

$$ \mathbf{K}_{10}^b = \{C_1(\mathbf{x}_{\alpha} - \mathbf{x}_0)\}, \alpha \in \mathcal{B}_{\alpha}. $$

$$ \mathbf{K}_{10}^n = \{C_1(\mathbf{x}_{\alpha} - \mathbf{x}_0)\}, \alpha \in \bar{\mathcal{B}}_{\alpha}. $$

$$ \mathbf{K}_{120}^b = \{C_{12}(\mathbf{x}_{\beta} - \mathbf{x}_0)\}, \beta \in \mathcal{B}_{\beta}. $$

$$ \mathbf{K}_{120}^n = \{C_{12}(\mathbf{x}_{\beta} - \mathbf{x}_0)\}, \beta \in \bar{\mathcal{B}}_{\beta}. $$

$$ \mathbf{g}_{2b}^j = \{g^j(\mathbf{x}_{\beta})\}, \beta \in \mathcal{B}_{\beta}, \mathbf{g}_{2n}^j = \{g^j(\mathbf{x}_{\beta})\}, \beta \in \bar{\mathcal{B}}_{\beta}. $$

$$ \mathbf{f}_{1b}^i = \{f^i(\mathbf{x}_{\alpha})\}, \alpha \in \mathcal{B}_{\alpha}, \mathbf{f}_{1n}^i = \{f^i(\mathbf{x}_{\alpha})\}, \alpha \in \bar{\mathcal{B}}_{\alpha}. $$

$$ \mathbf{f}_{2b}^i = \{f^i(\mathbf{x}_{\beta})\}, \beta \in \mathcal{B}_{\beta}, \mathbf{f}_{2n}^i = \{f^i(\mathbf{x}_{\beta})\}, \beta \in \bar{\mathcal{B}}_{\beta}. $$

$$ f_0^i = f^i(\mathbf{x}_0) \text{ is the ith covariate } f^i \text{ at the target location } \mathbf{x}_0. $$

Because the covariance function is assumed to be symmetric, see again Beauchamp et al. [3], $\mathbf{K}_{21}^{\text{bb}}$, $\mathbf{K}_{21}^{\text{bn}}$, $\mathbf{K}_{21}^{\text{nb}}$ and $\mathbf{K}_{21}^{\text{nn}}$ are respectively the transposes of the matrices $\mathbf{K}_{12}^{\text{bb}}$, $\mathbf{K}_{12}^{\text{bn}}$, $\mathbf{K}_{12}^{\text{nn}}$ and $\mathbf{K}_{12}^{\text{nn}}$.

$\mathbf{1}_{11}^{\text{bb}}$ and $\mathbf{0}_{11}^{\text{bb}}$ are the identity and zero matrices with same size than $\mathbf{K}_{11}^{\text{bb}}$. The same applies for all the 16 possible combinations of covariance matrices between basic and non-basic weights of $Y$ and $Z$.
Expanding the system (18a,18b,18c,18d) and because, by definition, the non-basic $\lambda_\alpha$ and $\nu_\beta$ respectively equals $\overline{\lambda_\alpha}$ and $\overline{\nu_\beta}$ and the basic $\tau_\alpha$ and $\omega_\beta$ equals 0:

$$ K_{11}^{bb} \lambda_b + K_{12}^{bb} \nu_b + \mu_0 + \nu_i f_{1b}^i = K_{10}^b - K_{11}^{bn} \overline{\lambda_b} - K_{12}^{bn} \overline{\nu_b} \quad (19a) $$

$$ K_{11}^{nb} \lambda_b + K_{12}^{nb} \nu_b + \mu_0 + \nu_i f_{1n}^i + 1^{nn} \tau_n = K_{10}^n - K_{11}^{nn} \overline{\lambda_b} - K_{12}^{nn} \overline{\nu_b} \quad (19b) $$

$$ K_{21}^{bb} \lambda_b + K_{22}^{bb} \nu_b + \nu_0 + \mu_j g_{2b}^j + \nu_i f_{2b}^i = K_{20}^b - K_{21}^{bn} \overline{\lambda_b} - K_{22}^{bn} \overline{\nu_b} \quad (19c) $$

$$ K_{21}^{nb} \lambda_b + K_{22}^{nb} \nu_b + \nu_0 + \mu_j g_{2n}^j + \nu_i f_{2n}^i + 1^{nn} \omega_n = K_{20}^n - K_{21}^{nn} \overline{\lambda_b} - K_{22}^{nn} \overline{\nu_b} \quad (19d) $$

$$ 1^b \lambda_b = 1 - 1^n \overline{\lambda_n} \quad (19e) $$

$$ g_{2b}^j \lambda_b = -g_{2n}^j \overline{\nu_n} \quad (19f) $$

$$ 1^b \nu_b = -1^n \overline{\nu_n} \quad (19g) $$

$$ f_{1b}^i \lambda_b + f_{2b}^j \nu_b = f_0^i - f_{1n}^i \lambda_n - f_{2n}^j \overline{\nu_n} \quad (19h) $$

$$ \lambda_\alpha - \overline{\lambda_\alpha} \le 0 \quad \text{and} \quad \tau_\alpha^* \ge 0, \quad \alpha = 1, \dots, N \quad (19i) $$

$$ \nu_\beta - \overline{\nu_\beta} \le 0 \quad \text{and} \quad \omega_\beta^* \ge 0, \quad \beta = 1, \dots, M $$

The condition (18c) is implicitly satisfied at optimality by the definitions of the basic and non-basic $\lambda$, $\nu$, $\tau$ and $\omega$.

Rewriting conditions (19a), (19c), (19e), (19f), (19g), (19h) in a convenient matrix form (non-basic conditions (19b) and (19d) are not embedded in the linear system) leads to:

$$ \begin{bmatrix} K_{11}^{bb} & K_{12}^{bb} & 1 & 0 & 0 & f_{1b}^i \\ K_{21}^{bb} & K_{22}^{bb} & 0 & g_{2b}^j & 1 & f_{2b}^i \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & g_{2b}^j & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ f_{1b}^i & f_{2b}^i & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \lambda_b \\ \nu_b \\ \mu_0 \\ \mu_j \\ \nu_0 \\ \nu_i \end{bmatrix} = \begin{bmatrix} K_{10}^b - K_{11}^{bn} \overline{\lambda_b} - K_{12}^{bn} \overline{\nu_b} \\ K_{20}^n - K_{21}^{bn} \overline{\lambda_b} - K_{22}^{bn} \overline{\nu_b} \\ 1 - 1^n \overline{\lambda_n} \\ -g_{2n}^j \overline{\nu_n} \\ -1^n \overline{\nu_n} \\ f_0^i - f_{1n}^i \overline{\lambda_n} - f_{2n}^i \overline{\nu_n} \end{bmatrix} $$

The final estimation of $Y(\mathbf{x}_0)$ will be:

$$ Y^{\text{PCK}}(\mathbf{x}_0) = \sum_{\alpha \in \mathcal{B}_\alpha} \lambda_\alpha Y(\mathbf{x}_\alpha) + \sum_{\alpha \in \overline{\mathcal{B}}_\alpha} \overline{\lambda}_\alpha Y(\mathbf{x}_\alpha) \\ + \sum_{\beta \in \mathcal{B}_\beta} \nu_\beta Z(\mathbf{x}_\beta) + \sum_{\beta \in \overline{\mathcal{B}}_\beta} \overline{\nu}_\beta Z(\mathbf{x}_\beta) \quad (20) $$

with $\mathcal{B}_\alpha$ and $\mathcal{B}_\beta$ the set of index for the basic weights of $Y(\mathbf{x}_\alpha)$ and $Z(\mathbf{x}_\beta)$, and $\overline{\mathcal{B}}_\alpha$ and $\overline{\mathcal{B}}_\beta$ their complementary, see Eq.(16) and (17).

The related cokriging variance of $Y^{\text{PCK}}(\mathbf{x}_0)$, Eq. (20), is (see the detailed calculations in AppendixA):

$$ \begin{aligned} \text{Var}[Y(\mathbf{x}_0) - Y^{\text{PCK}}(\mathbf{x}_0)] &= \sigma_{\text{PCK}}^2(\mathbf{x}_0) \\ &= C^1(\mathbf{0}) - \sum_{\alpha \in \mathcal{B}_\alpha} \lambda_\alpha C_{\alpha 0}^1 - \sum_{\alpha \in \overline{\mathcal{B}}_\alpha} \lambda_\alpha [C_{\alpha 0}^1 - \tau_\alpha] \\ &\quad - \sum_{\beta \in \mathcal{B}_\beta} \nu_\beta C_{\beta 0}^{12} - \sum_{\beta \in \overline{\mathcal{B}}_\beta} \nu_\beta [C_{\beta 0}^{12} - \omega_\beta] - \mu_0 - \sum_i \nu_i f_i^i(\mathbf{x}_0) \end{aligned} $$

where $\{\tau_\alpha\}$ and $\{\omega_\beta\}$ are directly given by Eq.(19b) and (19d).

## 2.4. Strategy

The basic and non-basic $\lambda_\alpha$ satisfy the condition (7), which requires to know the basic and non-basic $\nu_\beta$ first, satisfying the condition (10). The first step thus consists in solving the cokriging system looking at each step of the algorithm which $\nu_\beta$ is the most non-basic, see section 2.6 and following Eq. (21), until there is no non-basic $\nu_\beta$ anymore.

$$ \begin{bmatrix} K_{11}^{bb} & K_{12}^{bb} & 1 & 0 & 0 & f_{1b}^i \\ K_{21}^{bb} & K_{22}^{bb} & 0 & g_{2b}^j & 1 & f_{2b}^i \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & g_{2b}^j & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ f_{1b}^i & f_{2b}^i & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \lambda_b \\ \nu_b \\ \mu_0 \\ \mu_j \\ \nu_0 \\ \nu_i \end{bmatrix} = \begin{bmatrix} K_{10}^b - K_{12}^{bn} \overline{\nu_b} \\ K_{20}^n - K_{22}^{bn} \overline{\nu_b} \\ 1 \\ -g_{2n}^j \overline{\nu_n} \\ -1^n \overline{\nu_n} \\ f_0^i - f_{2n}^i \overline{\nu_n} \end{bmatrix} \quad (21) $$

Then, the $\nu_\beta$ are considered as known and fixed. The cok-
<figure><img src="image_3.png" /><figcaption>**Figure 1:** algorithm of the penalized cokriging (PCK) algorithm</figcaption></figure>

riging system is after that modified to satisfy the condition (7). The terms related to the second variable Z are moved to the right-hand side of the linear system, so that the unbiasedness conditions and optimality are still satisfied:

$$ \begin{bmatrix} \mathbf{K}_{11}^{\mathrm{bb}} & 1 & \mathbf{f}_{1\mathrm{b}}^{i} \\ 1 & 0 & 0 \\ \mathbf{f}_{1\mathrm{b}}^{i} & 0 & 0 \end{bmatrix} \begin{bmatrix} \mathbf{\lambda}_{\mathrm{b}} \\ \mu_{0} \\ \mathbf{v}_{i} \end{bmatrix} = \begin{bmatrix} \mathbf{K}_{10}^{\mathrm{b}} - \mathbf{K}_{11}^{\mathrm{bn}} \bar{\lambda}_{n} - \mathbf{K}_{12}^{\mathrm{bn}} \bar{\nu}_{n} - \mathbf{K}_{12}^{\mathrm{bb}} \bar{\nu}_{b} \\ 1 - \mathbf{1}^{n} \bar{\lambda}_{n} \\ f_{0}^{i} - \mathbf{f}_{1n}^{i} \bar{\lambda}_{n} - \mathbf{f}_{2n}^{i} \bar{\nu}_{n} - \mathbf{f}_{2b}^{i} \bar{\nu}_{b} \end{bmatrix} \quad (22) $$

The algorithm may still reduce the set of basic weights too far and the estimation will just be a heuristic produced by the additional constraints made on the cokriging system. Because two unbiasedness conditions are still appearing in the system (22), two basic weights are at least required.

## 2.5. Computational details

Given that A is symmetric, $A^{-1}$ exists and the inverse of the matrix A extended by an additional row and column is:

$$ \begin{bmatrix} \mathbf{A} & \mathbf{a} \\ \mathbf{a}' & \alpha \end{bmatrix}^{-1} = \begin{bmatrix} \mathbf{C} & \mathbf{c} \\ \mathbf{c}' & -\beta \end{bmatrix} $$

### subprocedures

isBasic : test if all the weights are basic

findMnb : find the most non-basic weight

updateCKDE: update the cokriging linear matrix system

updateKDE : update the kriging linear matrix system

estimPM<sub>2.5</sub> : compute the PM<sub>2.5</sub> penalized cokriging

### variables

$\bar{\lambda}, \bar{\nu}$ : weights of $Y(\mathbf{x}_\alpha)$ and $Z(\mathbf{x}_\beta)$ in $Y^{\mathrm{PCK}}(\mathbf{x}_0)$

$\bar{\lambda}$ : new constraints on $\lambda$

$\bar{\nu}$ : new constraints on $\nu$

$\lambda_b, \lambda_n$ : basic and non-basic $\lambda$

$\nu_b, \nu_n$ : basic and non-basic $\nu$

$\mathbf{K}$

$\mathbf{A}$

$\mathbf{W}$ : matrices and vectors to update

$\mathbf{K}_0$

$\mathbf{F}_0$

$Y^{\mathrm{PCK}}(\mathbf{x}_0)$: PM<sub>2.5</sub> estimation made by the algorithm (PCK)

$Z^{\mathrm{CK}}(\mathbf{x}_0)$: PM<sub>10</sub> estimation made by the algorithm (CK)

$$ \begin{bmatrix} \mathbf{K} & \mathbf{A} \\ \mathbf{A}' & \mathbf{0} \end{bmatrix}, \begin{bmatrix} \mathbf{W} \\ \mathbf{M} \end{bmatrix}, \begin{bmatrix} \mathbf{K}_0 \\ \mathbf{F}_0 \end{bmatrix} $$

Thus, $A^{-1}$ can be efficiently computed: $A^{-1} = C+(1/\beta)cc'$.

When classifying a weight $\lambda_\alpha$ or $\nu_\beta$ as non-basic, the row and the column to remove are generally not the last ones. The tool to compute the inverse of a matrix A after a permutation of the row and column j to the far right and bottom of matrix A is also given in [2]:

$$ \begin{bmatrix} a_{11} & \cdots & a_{1j} & \cdots & a_{1n} \\ \vdots & & \vdots & & \vdots \\ a_{j1} & \cdots & a_{jj} & \cdots & a_{jn} \\ \vdots & & \vdots & & \vdots \\ a_{n1} & \cdots & a_{nj} & \cdots & a_{nn} \end{bmatrix}^{-1} = \begin{bmatrix} c_{11} & \cdots & c_{1j} & \cdots & c_{1n} \\ \vdots & & \vdots & & \vdots \\ c_{j1} & \cdots & c_{jj} & \cdots & c_{jn} \\ \vdots & & \vdots & & \vdots \\ c_{n1} & \cdots & c_{nj} & \cdots & c_{nn} \end{bmatrix} \\ \Rightarrow \begin{bmatrix} a_{11} & \cdots & a_{1n} & a_{1j} \\ \vdots & & \vdots & \vdots \\ a_{n1} & \cdots & a_{nn} & a_{nj} \\ \vdots & & \vdots & \vdots \\ a_{j1} & \cdots & a_{jn} & a_{jj} \end{bmatrix}^{-1} = \begin{bmatrix} c_{11} & \cdots & c_{1n} & c_{1j} \\ \vdots & & \vdots & \vdots \\ c_{n1} & \cdots & c_{nn} & c_{nj} \\ \vdots & & \vdots & \vdots \\ c_{j1} & \cdots & c_{jn} & c_{jj} \end{bmatrix} \end{array} $$

Instead of solving the new linear system at each iteration of the algorithm, it is possible with the previously defined
tools to quickly update the inverse of the matrix and simply multiply it with the update of the second member. An even faster solution is to use some simple properties of linear algebra. Let us note $[y \ z]'$ the solution of the linear system at step K of the algorithm. Thus,

$$ \begin{bmatrix} \mathbf{A} & \mathbf{a} \\ \mathbf{a}' & \alpha \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{K} \\ k \end{bmatrix} = \begin{bmatrix} \mathbf{C} & \mathbf{c} \\ \mathbf{c}' & -\beta \end{bmatrix} \begin{bmatrix} \mathbf{K} \\ k \end{bmatrix} = \begin{bmatrix} \mathbf{y} \\ \mathbf{z} \end{bmatrix} $$

At step K + 1, a new observation is defined as non-basic. The new basic weights are obtained by multiplying the matrix $\mathbf{A}^{-1}$ with the column vector $\mathbf{K} - \mathbf{a}\bar{z}$, where $\bar{z}$ is the non-basic weights vector related to z, i.e.:

$$ \begin{aligned} \mathbf{A}^{-1}(\mathbf{K} - \mathbf{a}\bar{z}) &= \mathbf{C} + \left(\frac{1}{\beta}\right)\mathbf{c}\mathbf{c}'(\mathbf{K} - \mathbf{a}\bar{z}) \\ &= \mathbf{C}\mathbf{K} + \frac{1}{\beta}\mathbf{c}\mathbf{c}'\mathbf{K} - \mathbf{C}\mathbf{a}\bar{z} - \frac{1}{\beta}\mathbf{c}\mathbf{c}'\mathbf{a}\bar{z} \\ &= \mathbf{y} + \frac{z}{\beta}\mathbf{c} - \mathbf{C}\mathbf{a}\bar{z} - \frac{1}{\beta}\mathbf{c}\mathbf{c}'\mathbf{a}\bar{z} \end{aligned} $$

By noticing that:

$$ \begin{bmatrix} \mathbf{C} & \mathbf{c} \\ \mathbf{c}' & -\beta \end{bmatrix} \begin{bmatrix} \mathbf{A} & \mathbf{a} \\ \mathbf{a}' & \alpha \end{bmatrix}^{-1} = \mathbf{I} $$

implies that:

$$ \begin{bmatrix} \mathbf{C} & \mathbf{c} \\ \mathbf{c}' & -\beta \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \alpha \end{bmatrix} = 0, \text{ i.e. } \mathbf{C}\mathbf{a}\bar{z} = -\alpha\mathbf{c}\bar{z} $$

$$ \begin{bmatrix} \mathbf{c}' & -\beta \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \alpha \end{bmatrix} = 1, \text{ i.e. } \mathbf{d}'\mathbf{a}\bar{z} = (1 + \beta\alpha)\bar{z} $$

Finally,

$$ \begin{aligned} \mathbf{A}^{-1}(\mathbf{K} - \mathbf{a}\bar{z}) &= (\mathbf{C} + (\frac{1}{\beta})\mathbf{c}\mathbf{c}')(\mathbf{K} - \mathbf{a}\bar{z}) \\ &= \mathbf{y} + \mathbf{c} \left( \frac{z - \bar{z}(1 + \alpha(\beta + 1))}{\beta} \right) \end{aligned} $$

The computational cost is thus considerably reduced. Instead of multiplying a squared matrix of rank $N(B) = \#\{B_\alpha \cup B_\beta\}$, i.e. the number of basic weights, with a column vector of size $N(B)$, that is $N(B)^2$ multiplications and $[N(B) - 1]^2$ additions, the numbers of operations required to update the basic weights is now only $N(B)+3$ additions and $N(B)+3$ multiplications.

## 2.6. Remarks

A) The computational time of the method presented here may remain quite restricting when using a unique neighborhood because cokriging matrices have a size much bigger than kriging matrices. Thus, despite the previous tools introduced in Sect. 2.5 to avoid the solving of the linear system at each step of the algorithm, the update of the second member, the matrix cokriging and the weights may still be too costly for an operational algorithm. A solution to

speed up the estimation can be to consider a moving neighborhood and is fully detailed in Appendix B.

B) At each step of the algorithm, if there are still some weights that do not satisfy the constraints (7), the "most" non-basic weights $\lambda_\alpha$ have to be removed from the set of basic weights. The most non-basic index $\alpha$ can be defined by:

Eq. (24), as the weight that satisfies the least the constraints (i.e that maximizes their differences):

$$ \underset{\alpha}{\operatorname{argmax}} \{\lambda_{\alpha} - \bar{\lambda}_{\alpha} > 0\} \qquad (24) $$

or as Eq. (25), as the weight non satisfying the constraints (7) and maximizing the distance with the target location $\mathbf{x}_0$ where an estimation has to be made:

$$ \underset{\alpha}{\operatorname{argmax}} \{ \|\mathbf{x}_{\alpha} - \mathbf{x}_0\|, \lambda_{\alpha} - \bar{\lambda}_{\alpha} > 0 \} \qquad (25) $$

This second option allows the weights related to the observations in a close neighborhood of $\mathbf{x}_0$ to be truly estimated and not set to their non-basic related values when the algorithm successively reduces the set of basic weights. The same applies to find the most non-basic weights $v_\beta$.

## 3. Results

In this Section, we proposed an application based on the dataset used in Beauchamp et al. [3]. The specific additive modelling introduced in the above-mentioned work is used to compare both mappings and performance of the personalized cokriging (PCK) algorithm with the usual way of solving the cokriging system (CK). The background PM$_{10}$ and background PM$_{2.5}$ observational dataset is then identical to Beauchamp et al. [3], but extended from France to Europe, see Table 1, because the discontinuities discussed in Appendix B are more frequently seen when the domain is larger.

<table><caption>**Table 1:** Number of background monitoring sites (PM$_{10}$ and PM$_{2.5}$)</caption><thead><tr><th></th><th>Rural</th><th>Suburban</th><th>Urban</th></tr></thead><tbody><tr><td>PM<sub>10</sub></td><td>190</td><td>183</td><td>445</td></tr><tr><td>PM<sub>2.5</sub></td><td>63</td><td>51</td><td>188</td></tr></tbody></table>

The CHIMERE PM$_{10}$ and PM$_{2.5}$ simulations are used as covariates for $m_Z(\mathbf{x})$ and $m_Y(\mathbf{x})$. The model covers the AWM European domain of simulation ($[-15°W, 35°W; 35°N, 70°N]$), with a coarse resolution of 0.5°. The maps are interpolated by cokriging on a regular grid with the same resolution used by CHIMERE for solving the chemical and physical processes.

The time period covers the first quarter of 2015 in which the PM concentrations were particularly high with a long episode of pollution occuring in March 2015.
Regarding the practical implementation, an original C++
program has been developed to run the penalized cokriging
version. The armadillo C++ library [12] is involved for the
numerical analysis related to kriging. An interface with R
software [29] is used inside the C++ program towards the
R library RcppArmadillo [17] so that all the descriptive
geostatistical part, e.g. the variogram computation and fit-
ting, are done by RGeostats through this interface. The full
code is available via Zenodo (10.5281/zenodo.7756425). It
is ready-to-use after installation of the appropriate libraries
and specifications of the input data. The appropriate shape
for the inputs is given as comment lines along the program.

**3.1. Mapping**

On the 10^th^ of March, Fig. 2 (a) to Fig. 2 (j) respec-
tively show the daily observations, the corresponding scat-
terplot between the collocated PM<sub>2.5</sub> and PM<sub>10</sub> sites, the
CHIMERE simulations for the two pollutants and the krig-
ing maps with their related standard deviations for PM<sub>10</sub>
and PM<sub>2.5</sub> with the classic algorithm (CK) and its penalized
version (PCK). A zoomed-in window is given on Hungary,
Serbia, Ukraine and Slovakia to focus on this area where
Y<sup>CK</sup>(x<sub>0</sub>) > Z<sup>CK</sup>(x<sub>0</sub>). The estimations Y<sup>PCK</sup>(x<sub>0</sub>) are suc-
cessfully less than Z<sup>CK</sup>(x<sub>0</sub>). Looking at the daily-averaged
observations in this area, that are only two in Hungary, the
algorithm (PCK) seems to extend the representativeness ar-
eas of these two monitoring sites over the whole area where
the inconsistency appears. PM<sub>10</sub> observations are in the
range 30-40 µg.m<sup>-3</sup> while PM<sub>2.5</sub> cokriging estimates val-
ues are greater than 40 µg.m<sup>-3</sup> over the area. Adding the
penalization enables to decrease PM<sub>2.5</sub> estimations down
to 30-40 µg.m<sup>-3</sup> which is not unrealistic when regarding
how the values are spatially distributed. Though, an ad-
ditional assessment of the true representativeness areas of
these sites would be necessary to conclude. The same type
of results is shown on the 16^th^ of March in Fig. 3 (a) to
Fig. 3 (j), that is a typical example of far-off extrapola-
tions in Northern Europe where physical inconsistencies
can also occur. Once again, the algorithm (PCK) is suc-
cessful and enables to decrease the prior values produced
by the algorithm (CK) from 30 µg.m<sup>-3</sup> to less than 10
µg.m<sup>-3</sup> in its (PCK) version. It is to note that standard
deviation of (PCK) errors were expected to be greater than
those produced without the penalization but in most cases,
as supported by the two examples given, they are in the
same range of values for both algorithms.

More generally, the differences are often not significant from
a mapping perspective. It was expected since the usual cok-
riging algorithm does not generate so many PM2.5 estima-
tions that are greater than their corresponding PM10 esti-
mations. On some days however, the differences between
the cokriging and its penalized version are significant, es-
pecially because the algorithm (PCK) has the direct con-
sequence to fix the (non-basic) weights of a some obser-
vations in the estimation process. As a consequence, even

if all the available observations are first used as input data,
the progressive decrease of the set of basic weights leads to
a map that seems to be built with a moving neighborhood.
This one has very specific features strongly depending on
the strategy used to define the most non-basic weight, see
Sect. 2.6:

1) if it is defined as the farthest non-basic weight from the
target point x₀, see Eq.(25), then all the observations with
basic weights are included in some distance-based neigh-
borhood. If this distance is large enough, which is generally
the case, the approach proposed in Appendix B makes van-
ish the discontinuities because the farthest observations, al-
though with fixed (non-basic) weights, are noisy when es-
timating the PM₁₀ concentration Z(x₀).

2) if it is defined as the most non-basic weight, i.e. the
weight with the largest positive deviation to its related non-
basic version, see Eq.(24), then the set of basic weights
fail to comply with some distance-based neighborhood. As
a consequence, the estimation process still may generate
strong discontinuities. It would clearly be the best option
to keep the largest set of basic weights, but for mapping
concerns, the way of dealing with these discontinuities ap-
pears problematic and Appendix B will not really help.

**3.2. Cross-validation**

Figure 5 shows the results of a (leave-one-out) cross-validation procedure carried out every day of the first quarter of 2015. Figure 5a presents the scatterplot of the PM<sub>2.5</sub> observations (X-axis) and the PM<sub>2.5</sub> estimations (Y-axis) obtained by the cokriging algorithm (CK) of Beauchamp et al. [3] and its penalized version (PCK) detailed in Sect. 2. In the latter, the definition (25) of what should be the most non-basic weight is used. Over the period, only the observations Y(x<sub>α</sub>) with inconsistent cross-validation values Y<sup>CK</sup>(x<sub>α</sub>) > Z<sup>CK</sup>(x<sub>α</sub>) are displayed and used to compute the correlation and RMSE. Let precise that observations sites with inconsistent cross-validation values have a significant number of occurences, see Figure 4 focusing on the French domain and same evaluation period, with related inconsistencies on classic cokriging mapping (without additive model) also noticeable, see the background mapping on the above-mentioned Figure.
<figure><img src="image_4.png" /><figcaption>**Figure 2:** From top to bottom, PM₁₀ and PM₂.₅ available observations and corresponding scatterplots, PM₁₀ and PM₂.₅ CHIMERE outputs, PM₁₀ cokriging and standard deviation as implemented in [3], the same for PM₂.₅ and last, the proposed PM₂.₅ penalized cokriging and standard deviation as proposed in this paper (2015, March 10). A focus is provided on a small box region over Hungary, Serbia, Ukraine and Slovakia with inconsistent cokriging estimations</figcaption></figure>
An iterative optimization scheme to accommodate inequality constraints in air quality geostatistical estimation of multivar. PM

<figure><img src="image_5.png" /><figcaption>**Figure 3:** From top to bottom, PM₁₀ and PM₂.₅ available observations and corresponding scatterplots, PM₁₀ and PM₂.₅ CHIMERE outputs, PM₁₀ cokriging and standard deviation as implemented in [3], the same for PM₂.₅ and last, the proposed PM₂.₅ penalized cokriging and standard deviation as proposed in this paper (2015, March 16). A focus is provided on a small box region o the western coastal area of Norway with inconsistent cokriging estimations</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Cross-validation procedure over the first quarter of 2015</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Frequencies of PM<sub>2.5</sub>/PM<sub>10</sub> ratio greater than 1 for both observation sites (cross-validation) and cokriging estimation maps</figcaption></figure>

In the case of consistent cokriging-based estimations, the scores of the two methods are similar.

$$
\begin{align*}
\text{1) correlation} \quad &: \frac{\text{Cov}[(Z(\mathbf{x}_{\alpha}), Z^{\text{CK}}(\mathbf{x}_{\alpha}))]}{\sigma[Z(\mathbf{x}_{\alpha})]\sigma[Z^{\text{CK}}(\mathbf{x}_{\alpha})]} \\
\text{2) (RMSE)} \quad &: \sqrt{\frac{1}{n} \sum_{\alpha=1}^{n} [Z(\mathbf{x}_{\alpha}) - Z^{\text{CK}}(\mathbf{x}_{\alpha})]^2}
\end{align*}
$$

The preliminary findings from Sect. 3.1 are now confirmed
by the validation procedure, which clearly shows that:

$$
Y^{\text{PCK}}(\mathbf{x}_{\alpha}) \leq Z^{\text{CK}}(\mathbf{x}_{\alpha}) \leq Y^{\text{CK}}(\mathbf{x}_{\alpha}),
$$

meaning that the iterative algorithm is successful. Regard-
ing the scores, if the correlation is a bit less with the new
algorithm, its RMSE is better. As a consequence, not only
the estimations satisfy the inequality constraint but they are
also consistent with the observations.

From the cross-validation procedure described previously,
we respectively store the values of #{B_α} and #{B_β}, the
number of basic weights for the main and secondary vari-
able. A high number of basic weights indicates that only
a few number of weights are fixed in the iterative process
so that the PM2.5 estimations satisfy the physical inequal-
ity. On the contrary, extremely low values for the quanti-
ties means that most of the observations available are used
with their non-basic weighting in the estimation, leading to
a simple heuristic. Figure 5b displays the bivariate distri-
bution of the basic weights. As we can see, only a few sam-
ples display low values for both #{B_α} and #{B_β}. Even
when it is the case, and because we use definition (25) for
the most non-basic weights identification, the nearest ob-
servations are always used to compute the PM2.5 optimal
interpolation, which explains why the estimations remain
consistent. In terms of computation cost, let precise that
because the number of modified weights are small, the pe-
nalized version cost is very close to the original cokriging
algorithm. When the inequality constraint is not satisfied
and the iterative penalized scheme is involved, less than
10 iterations were generally involved in our datasets (i.e.
less than 10 observational weights are non-basic) and no
estimation location ended with a simple heuristic interpo-
lation (when the number of iteration reaches the number of
observations, thus excluding most of the available informa-
tion). This specific cases happened when PM2.5 and PM10
observations are both close in terms of values and spatial
locations.

4. Conclusions

The problem of consistency in cokriging arises when deal-
ing with quantities that involve inequality constraints. In
a study by Beauchamp et al. [3], cokriging was employed to enhance the estimation of PM<sub>2.5</sub> by incorporating PM<sub>10</sub> observations using a specific additive modeling approach. The local means were derived by exploiting the physical relationship between the two variables. Although this cokriging approach led to improved estimations, it failed to ensure that the resulting estimations adhered to the inequality constraint. From a mathematical standpoint, it is possible to verify this by conducting conditional simulations at a target location x<sub>0</sub>. Such simulations reveal that the intersection of the PM<sub>2.5</sub> simulated distribution with the PM<sub>10</sub> simulated distribution is never empty. Thus, even when the average PM<sub>2.5</sub> simulation is higher than the average PM<sub>10</sub> simulation, the consistency of the estimation is preserved.

To address the issue of inconsistencies in PM<sub>2.5</sub> concentrations, a new algorithm is proposed, in which additional constraints are introduced on the cokriging weights. They allow for successive iterations to solve the cokriging system in terms of basic component that satisfies these new constraints. The algorithm also includes several computational details that ensure a reasonable computational cost, making it practical for operational contexts.

Although the new estimator is performing well in most cases, there are still some limitations that need to be considered. One such limitation is the algorithm sensitivity to observation noise, as inaccuracies in the observed data can lead to significant errors in the final mapping. This can be especially problematic in areas with a high degree of spatial variability, where the data may be sparse or irregularly distributed. In such cases, the successive iterations in the penalization procedure may reduce the set of basic weights too far and thus build a simple heuristic that poorly estimates the true PM<sub>2.5</sub> concentration. This can result in a less accurate final mapping and may require additional adjustments to improve the algorithm's performance. Other limitations of the proposed algorithm include its reliance on assumptions about the underlying spatial structure of the data, as well as its potential limitations in handling non-Gaussian and non-stationary data, which could be crucial for communication and decision-making.

Finally, the use of our algorithm may apply to similar problems: it would remain identical and valid. First example would be the estimation of PM non-volatile fraction, see e.g. [5]. Only the penalization $\bar{\lambda}$ and $\bar{v}$ introduced on the cokriging weights shall be adapted.

### Acknowledgements

This study was partially funded by the French Ministry in charge of the environment, when the authors were members of the French Central Laboratory for Air Quality Monitoring (LCSQA).

### Author contribution statement

**Maxime Beauchamp:** Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper. **Bertrand Bessagnet:** Analyzed and interpreted the data; Wrote the paper.

### Data availability statement

The software developed is now available through a DOI produced via Zenodo. The dataset comes from a public database (CAMS Copernicus) but since this work was made 4-5 years it was not possible to fully retrieve the exact dataset used at this time. We prefer to make the code public with associated data formatting so that any potential user can adapt and test the code on new datasets.

### Declaration of interest's statement

The authors declare no conflict of interest.
References

[1] P. Abrahamsen and F. E. Benth. Kriging with inequality constraints. *Mathematical Geology*, 33(6):719–744, 2001. ISSN 1573-8868. doi: 10.1023/A:1011078716252. URL http://dx.doi.org/10.1023/A:1011078716252.

[2] R.J. Barnes and T.B. Johnson. Positive kriging. In G. Verly, M. David, A.G. Journel, and A. Marechal, editors, *Geostatistics for Natural Resources Characterization*, pages 231–244. Springer Netherlands, 1984. ISBN 978-94-010-8157-3. doi: 10.1007/978-94-009-3699-7_14. URL http://dx.doi.org/10.1007/978-94-009-3699-7_14.

[3] M. Beauchamp, B. Bessagnet, C. de Fouquet, L. Malherbe, F. Meleux, and A. Ung. An additive geostatistical model for mixing total and partial PM₁₀ observations with CHIMERE rCTM. *Atmospheric Environment*, 189:61–79, 2018. ISSN 1352-2310. doi: https://doi.org/10.1016/j.atmosenv.2018.06.035. URL http://www.sciencedirect.com/science/article/pii/S1352231018304266.

[4] B. Bessagnet, E. Debry, and L. Malherbe. Bilan des mesures de PM₁₀ et PM₂.₅ ajustées et évaluation des outils de modélisation. Technical Report DRC-10-111613-05227A, LCSQA, 2010. URL www.lcsqa.org.

[5] Bertrand Bessagnet, Nadine Allemand, Jean-Philippe Putaud, Florian Couvidat, Jean-Marc André, David Simpson, Enrico Pisoni, Benjamin N Murphy, and Philippe Thunis. Emissions of carbonaceous particulate matter and ultrafine particles from Vehicles-A scientific review in a Cross-Cutting context of air pollution and climate change. *Appl Sci (Basel)*, 12(7):1–52, April 2022.

[6] S. Boyd and L. Vandenberghe. *Convex optimization*. Cambridge University Press, Cambridge, 2004. ISBN 0-521-83378-7.

[7] P. J. Brown, N.D. Le, and J.V. Zidek. Multivariate spatial interpolation and exposure to air pollutants. *Canadian Journal of Statistics*, 22(4):489–509, 1994. ISSN 1708-945X. doi: 10.2307/3315406. URL http://dx.doi.org/10.2307/3315406.

[8] M. Cameletti, F. Lindgren, D. Simpson, and H. Rue. Spatiotemporal modeling of particulate matter concentration through the spde approach. *AStA Advances in Statistical Analysis*, 97(2):109–131, 2012. ISSN 1863-818X. doi: 10.1007/s10182-012-0196-3. URL http://dx.doi.org/10.1007/s10182-012-0196-3.

[9] J-P. Chilès, S. Séguret, and P-M. Riboud. Geostatistical analysis of validation data of an air pollution simulator. In *8th International geostatistics Congress*, page on CD, Santiago, Chile, December 2008. URL https://hal-mines-paristech.archives-ouvertes.fr/hal-00776934.

[10] J.P. Chiles and P. Delfiner. *Geostatistics: modeling spatial uncertainty*. Wiley, New-York, second edition, 2012.

[11] CIESIN. Global annual average PM₂.₅ grids from modis and misr aerosol optical depth (aod), 20160510 2013. URL http://dx.doi.org/10.7927/H4H41PB4.

[12] S. Conrad. Armadillo: an open source C++ linear algebra library for fast prototyping and computationally intensive experiments. Technical report, NICTA, 2010.

[13] N. Cressie and C. K. Wikl. The variance-based cross-variogram you can add apples and oranges. *Mathematical Geology*, 30:789–800, 1998.

[14] C. V. Deutsch. Correcting for negative weights in ordinary kriging. *Comput. Geosci.*, 22(7):765–773, August 1996. ISSN 0098-3004. doi: 10.1016/0098-3004(96)00005-2. URL http://dx.doi.org/10.1016/0098-3004(96)00005-2.

[15] Directive 2008/50/EC. *Directive 2008/50/EC of the European Parliament and of the Council of 21 May 2008 on ambient air quality and cleaner air for Europe*, 2008.

[16] O. Dubrule and C. Kostov. An interpolation method taking into account inequality constraints: I. methodology. *Mathematical Geology*, 18(1):33–51, 1986.

[17] D. Eddelbuettel and C. Sanderson. Rcpparmadillo: Accelerating R with high-performance C++ linear algebra. *Computational Statistics and Data Analysis*, 71:1054–1063, March 2014. URL http://dx.doi.org/10.1016/j.csda.2013.02.005.

[18] A.E. Gelfand, A.M. Schmidt, S. Banerjee, and C.F. Sirmans. Non-stationary multivariate process modeling through spatially varying coregionalization. *Test*, 13(2):263–312, Dec 2004. ISSN 1863-8260. doi: 10.1007/BF02595775. URL https://doi.org/10.1007/BF02595775.

[19] L. E. Gerharz and E. Pebesma. Using geostatistical simulation to disaggregate air quality model results for individual exposure estimation on gps tracks. *Stochastic Environmental Research and Risk Assessment*, 27(1):223–234, 2013. ISSN 1436-3259. doi: 10.1007/s00477-012-0578-9. URL http://dx.doi.org/10.1007/s00477-012-0578-9.

[20] G. B. M. Heuvelink, E. Pebesma, and B. Gräler. *Encyclopedia of GIS*, chapter Space-Time Geostatistics, pages 1–7. Springer International Publishing, Cham, 2015. ISBN 978-3-319-23519-6. doi: 10.1007/978-3-319-23519-6_1647-1. URL http://dx.doi.org/10.1007/978-3-319-23519-6_1647-1.

[21] E. H. Isaaks and R. M. Srivastava. *An Introduction to Applied Geostatistics*. Oxford University Press, New-York, 1989. ISBN 0-19-505012-6.

[22] K. Krivoruchko and A. Gribov. *geoENV IV — Geostatistics for Environmental Applications: Proceedings of the Fourth European Conference on Geostatistics for Environmental Applications held in Barcelona, Spain, November 27–29, 2002*, chapter Geostatistical Interpolation and Simulation in the Presence of Barriers, pages 331–342. Springer Netherlands, Dordrecht, 2004. ISBN 978-1-4020-2115-2. doi: 10.1007/1-4020-2115-1_28. URL http://dx.doi.org/10.1007/1-4020-2115-1_28.

[23] V. Langlais. Krigeage sous contraintes d'inégalités en voisinage glissant. Technical Report Note interne, Centre de Géostatistique, ENSMP, Fontainebleau, 1988.

[24] S. Mailler, L. Menut, D. Khvorostyanov, M. Valari, F. Couvidat, G. Siour, S. Turquety, R. Briant, P. Tuccella, B. Bessagnet, A. Colette, L. Létinois, K. Markakis, and F. Meleux. Chimere-2017: from urban to hemispheric chemistry-transport modeling. *Geoscientific Model Development*, 10(6):2397–2423, 2017. doi: 10.5194/gmd-10-2397-2017. URL https://www.geosci-model-dev.net/10/2397/2017/.

[25] S. Mailler, L. Menut, D. Khvorostyanov, M. Valari, F. Couvidat, G. Siour, S. Turquety, R. Briant, P. Tuccella, B. Bessagnet, A. Colette, L. Létinois, K. Markakis, and F. Meleux. Chimere-2017: from urban to hemispheric chemistry-transport modeling. *Geoscientific Model Development*, 10(6):2397–2423, 2017. doi: 10.5194/gmd-10-2397-2017. URL https://www.geosci-model-dev.net/10/2397/2017/.

[26] L. Malherbe and A. Ung. Travaux relatifs à la plate-forme nationale de modélisation PREV'AIR : Réalisation de cartes analysées d'Ozone. Technical Report DRC-10-103351-01139A, LCSQA, 2009. URL www.lcsqa.org.

[27] A. M. Michalak. A gibbs sampler for inequality-constrained geostatistical interpolation and inverse modeling. *Water Resources Research*, 44(9):n/a–n/a, 2008. ISSN 1944-7973. doi: 10.1029/2007WR006645. URL http://dx.doi.org/10.1029/2007WR006645. W09437.

[28] A. M. Michalak and P. K. Kitanidis. A method for enforcing parameter nonnegativity in bayesian inverse problems with an application to contaminant source identification. *Water Resources Research*, 39 (2):n/a–n/a, 2003. ISSN 1944-7973. doi: 10.1029/2002WR001480. URL http://dx.doi.org/10.1029/2002WR001480. 1033.

[29] R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2021. URL https://www.R-project.org/.

[30] M. C. Ribeiro, P. Pinho, C. Branquinho, E. Llop, and M. J. Pereira. Geostatistical uncertainty of assessing air quality using high-spatial-resolution lichen data: A health study in the urban area of sines, portugal. *Science of The Total Environment*, 562:740–750, 2016. ISSN 0048-9697. doi: http://dx.doi.org/10.1016/j.scitotenv.2016.04.081. URL http://www.sciencedirect.com/science/article/pii/
[31] J. Rivoirard and T. Romary. Continuity for kriging with moving neighborhood. *Mathematical Geosciences*, 43(4):469–481, 2011. ISSN 1874-8953. doi: 10.1007/s11004-011-9330-0. URL http://dx.doi.org/10.1007/s11004-011-9330-0.

[32] R. I. Rothenberg. *Linear Programming*. Oxford, North-Holland, New York, 1979.

[33] L. Rouil, C. Honoré, R. Vautard, M. Beekmann, B. Bessagnet, L. Malherbe, F. Meleux, A. Dufour, C. Elichegaray, J.-M. Flaud, L. Menut, D. Martin, A. Peuch, V.-H. Peuch, and N. Poisson. Prev’air: An operational forecasting and mapping system for air quality in europe. *Bulletin of the American Meteorological Society*, 90(1):73–84, 2009. doi: 10.1175/2008BAMS2390.1. URL https://doi.org/10.1175/2008BAMS2390.1.

[34] P.D. Sampson, M. Richards, A. A. Szpiro, S. Bergen, L. Shepard, T. V. Larson, and J. D. Kaufman. A regionalized national universal kriging model using partial least squares regression for estimating annual PM<sub>2.5</sub> concentrations in epidemiology. *Atmospheric Environment*, 75:383–392, 2013. ISSN 1352-2310. doi: http://dx.doi.org/10.1016/j.atmosenv.2013.04.015. URL http://www.sciencedirect.com/science/article/pii/S1352231013002604.

[35] A.M. Schmidt and A.E. Gelfand. A bayesian coregionalization approach for multivariate pollutant data. *Journal of Geophysical Research: Atmospheres*, 108(D24), 2003. ISSN 2156-2202. doi: 10.1029/2002JD002905. URL http://dx.doi.org/10.1029/2002JD002905. 8783.

[36] S. Soltani-Mohammadi and A. E. Tercan. Constrained multiple indicator kriging using sequential quadratic programming. *Computers & Geosciences*, 48:211–219, 2012. ISSN 0098-3004. doi: https://doi.org/10.1016/j.cageo.2012.01.003. URL http://www.sciencedirect.com/science/article/pii/S0098300412000052.

[37] F. Szidarovszky, E. Y. Baafi, and Y. C. Kim. Kriging without negative weights. *Mathematical Geology*, 19(6):549–559, Aug 1987. ISSN 1573-8868. doi: 10.1007/BF00896920. URL https://doi.org/10.1007/BF00896920.

[38] H. Wackernagel. *Multivariate Geostatistics. An Introduction with Applications*. Springer-Verlag Berlin Heidelberg, New York, 2003. doi: 10.1007/978-3-662-05294-5. ISBN 978-3-540-44142-7.

[39] Y. Yang, G. Christakos, W. Huang, C. Lin, P. Fu, and Y. Mei. Uncertainty assessment of PM<sub>2.5</sub> contamination mapping using spatiotemporal sequential indicator simulations and multi-temporal monitoring data. *Scientific Reports*, 6:24335 EP –, Apr 2016. URL http://dx.doi.org/10.1038/srep24335. Article.

# Appendices

## A. Variance of the penalized cokriging estimator

The related cokriging variance of $Y^{\text{PCK}}(\mathbf{x}_0)$ is:

$$ 
\begin{aligned}
& \text{Var}[Y(\mathbf{x}_0) - Y^{\text{PCK}}(\mathbf{x}_0)] = \sigma_{\text{PCK}}^2(\mathbf{x}_0) \\
& = C^Y(0) + \left\{ \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \alpha' \in \bar{B}_\alpha}} \lambda_\alpha \lambda_{\alpha'} C_{\alpha \alpha'}^1 + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \alpha' \in \bar{B}_\alpha}} \lambda_\alpha \lambda_{\alpha'} C_{\alpha \alpha'}^1 \right. \\
& \qquad \left. + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \alpha' \in \bar{B}_\alpha}} \lambda_\alpha \lambda_{\alpha'} C_{\alpha \alpha'}^1 + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \alpha' \in \bar{B}_\alpha}} \lambda_\alpha \lambda_{\alpha'} C_{\alpha \alpha'}^1 \right\} \\
& + \left\{ \sum_{\substack{\beta \in \bar{B}_\beta \\ \beta' \in \bar{B}_\beta}} \nu_\beta \nu_{\beta'} C_{\beta \beta'}^2 + \sum_{\substack{\beta \in \bar{B}_\beta \\ \beta' \in \bar{B}_\beta}} \nu_\beta \nu_{\beta'} C_{\beta \beta'}^2 \right. \\
& \qquad \left. + \sum_{\substack{\beta \in \bar{B}_\beta \\ \beta' \in \bar{B}_\beta}} \nu_\beta \nu_{\beta'} C_{\beta \beta'}^2 + \sum_{\substack{\beta \in \bar{B}_\beta \\ \beta' \in \bar{B}_\beta}} \nu_\beta \nu_{\beta'} C_{\beta \beta'}^2 \right\} \\
& - 2 \left\{ \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \alpha' \in \bar{B}_\alpha}} \lambda_\alpha C_{\alpha 0}^1 + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \alpha' \in \bar{B}_\alpha}} \lambda_\alpha C_{\alpha 0}^1 \right. \\
& \qquad \left. + \sum_{\substack{\beta \in \bar{B}_\beta \\ \beta' \in \bar{B}_\beta}} \nu_\beta C_{\beta 0}^1 + \sum_{\substack{\beta \in \bar{B}_\beta \\ \beta' \in \bar{B}_\beta}} \nu_\beta C_{\beta 0}^1 \right\} \\
& + 2 \left\{ \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \beta \in \bar{B}_\beta}} \lambda_\alpha \nu_\beta C_{\alpha \beta}^1 + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \beta \in \bar{B}_\beta}} \lambda_\alpha \nu_\beta C_{\alpha \beta}^1 \right. \\
& \qquad \left. + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \beta \in \bar{B}_\beta}} \lambda_\alpha \nu_\beta C_{\alpha \beta}^1 + \sum_{\substack{\alpha \in \bar{B}_\alpha \\ \beta \in \bar{B}_\beta}} \lambda_\alpha \nu_\beta C_{\alpha \beta}^1 \right\}
\end{aligned}
$$

By regrouping all the terms from Eq. (19a) to (19d), this
kriging variance is:

$$
\begin{align*}
\sigma_{\text{PCK}}^2(\mathbf{x}_0) \\
&= C^Y(\mathbf{0}) \\
&\quad - \sum_{\alpha \in \mathcal{B}_\alpha} \lambda_\alpha C_{\alpha 0}^1 - (\mu_0 + \sum_i v_i f^i(\mathbf{x}_\alpha)) \\
&\quad - \sum_{\alpha \in \mathcal{B}_\alpha} \lambda_\alpha [C_{\alpha 0}^1 - \tau_\alpha] - (\mu_0 + \sum_i v_i f^i(\mathbf{x}_\alpha)) \\
&\quad - \sum_{\beta \in \mathcal{B}_\beta} v_\beta C_{\beta 0}^{12} - (v_0 + \sum_i v_i f^i(\mathbf{x}_\beta) + \sum_j \mu_j g^j(\mathbf{x}_\beta)) \\
&\quad - \sum_{\beta \in \mathcal{B}_\beta} v_\beta [C_{\beta 0}^{12} - \omega_\beta] - (v_0 + \sum_i v_i f^i(\mathbf{x}_\beta) \\
&\quad + \sum_j \mu_j g^j(\mathbf{x}_\beta))
\end{align*}
$$

finally leading to:

$$
\begin{align*}
\sigma_{\text{PCK}}^2(\mathbf{x}_0) = C^1(\mathbf{0}) & - \sum_{\alpha \in \mathcal{B}_\alpha} \lambda_\alpha C_{\alpha 0}^1 - \sum_{\alpha \in \mathcal{B}_\alpha} \lambda_\alpha [C_{\alpha 0}^1 - \tau_\alpha] \\
& - \sum_{\beta \in \mathcal{B}_\beta} v_\beta C_{\beta 0}^{12} - \sum_{\beta \in \mathcal{B}_\beta} v_\beta [C_{\beta 0}^{12} - \omega_\beta] \\
& - \mu_0 - \sum_i v_i f^i(\mathbf{x}_0)
\end{align*}
$$

B. Continuous cokriging for moving neighborhood

In the usual cokriging framework, a moving neighborhood is more costly than a unique neighborhood because the kriging matrix has to be inverted for each target location x₀. In our penalized algorithm, because the system is updated at each step of the algorithm, the smaller this linear system is, the faster is the update. In addition, a moving neighborhood allows to refine at best the coefficients of the drift according to the local behavior of the covariates to the observations. Still, one issue of using a moving neighborhood is than it can create artificial discontinuities in the estimation when mapping the field on a (regular) grid.

In Rivoirard and Romary [31], the discontinuities caused by the moving neighbourhood are managed by considering the observations Z(x_β) spoiled by a noise ε_Z(x_β).

As a consequence, the kriging matrix is modified: to each term Cov[Z(x_β), Z(x_β')] is added Cov[ε_Z(x_β), ε_Z(x_β')]. The same applies for Cov[Y(x_α), Y(x_α')] and the kriging fashion is easily transposed for multivariate datasets. In particular, the cokriging becomes:

$$
\begin{align*}
Y^{\text{CK}}(\mathbf{x}_0) = & \sum_{\alpha} \lambda_{\alpha} \{Y(\mathbf{x}_{\alpha}) + \varepsilon_{Y}(\mathbf{x}_{\alpha})\} \\
& + \sum_{\beta} v_{\beta} \{Z(\mathbf{x}_{\beta}) + \varepsilon_{Z}(\mathbf{x}_{\beta})\}
\end{align*}
$$

When the initial covariance structure is non-continuous (with a nugget effect), ε_Z(x_β) and ε_Z(x_β') are considered independent when β ≠ β' and Cov[ε_Z(x_β), ε_Z(x_β')] = 0. The same applies for ε_Y. The variances Var[ε_Z(x_α)] and Var[ε_Y(x_β)] of the noises ε_Z and ε_Y are chosen by the user: they increases according to the distance to x_0 and are neglected or even set to 0 for the nearest data points:

$$
\begin{align*}
\operatorname{Var}[\varepsilon_Z(\mathbf{x}_\alpha)] &= C_Z(0) \left( \frac{\mathbf{h}_{\alpha 0} - r}{R - \mathbf{h}_{\alpha 0}} \right)^2 \\
\operatorname{Var}[\varepsilon_Y(\mathbf{x}_\beta)] &= C_Y(0) \left( \frac{\mathbf{h}_{\beta 0} - r}{R - \mathbf{h}_{\beta 0}} \right)^2
\end{align*}
$$

with h_{α0} = ||x_α - x_0||, h_{β0} = ||x_β - x_0||, R the radius of the moving neighborhood and r = p.R, with p < 1.

The lagrangian L is now defined as:

$$
\begin{align*}
L(\mathbf{x}_0) ={}& \operatorname{Var}[Y(\mathbf{x}_0) - Y^{\text{CK}}(\mathbf{x}_0)] + 2\mu_0 \left( \sum_{\alpha} \lambda_{\alpha} - 1 \right) \\
& + 2v_0 \sum_{\beta} v_{\beta} + 2\mu_j \left( \sum_{\alpha} \lambda_{\alpha} g^{j}(\mathbf{x}_{\alpha}) - g_0^j \right) \\
& + 2v_i \left( \sum_{\beta} v_{\beta} f^{i}(\mathbf{x}_{\beta}) - f_0^i \right) + \sum_{\alpha} \lambda_{\alpha}^2 n_{\alpha} + \sum_{\beta} v_{\beta}^2 n_{\beta}
\end{align*}
$$

where the weights λ_α and v_β are respectively penalized by
the quantities n_α = Var[ε_Y(x_β)] and n_β = Var[ε_Z(x_α)].

Equating the partial derivatives $\frac{\partial L}{\partial \lambda(\mathbf{x}_\alpha)}$ and $\frac{\partial L}{\partial v(\mathbf{x}_\beta)}$ to zero

leads to:

$$
\begin{align*}
& \lambda_{\alpha} n_{\alpha} + \sum_{\alpha'} \lambda_{\alpha'} C_{\alpha \alpha'}^{Y} - C_{\alpha 0}^{Y} + \sum_{\beta} v_{\beta} C_{\alpha \beta}^{12} \\
& \qquad + \mu_{0} + v_{i} f^{i}(\mathbf{x}_{\alpha}) + \mu_{j} g^{j}(\mathbf{x}_{\alpha}) &&= 0
\end{align*}
$$

and

$$
\begin{align*}
& v_{\beta} n_{\beta} + \sum_{\beta'} v_{\beta'} C_{\beta \beta'}^{2} - C_{\beta 0}^{12} + \sum_{\alpha} \lambda_{\alpha} C_{\alpha \beta}^{12} \\
& \qquad + v_{0} + v_{i} f^{i}(\mathbf{x}_{\beta}) &&= 0
\end{align*}
$$

and the cokriging matrix system is:

$$
(K + N)W + AM = K_0
$$

$$
A'W = F_0
$$

$$
\text{where } \mathbf{N} = \mathbf{In}' \text{ and } \mathbf{n} = \begin{bmatrix} \mathbf{n}_{\alpha} \\ \mathbf{n}_{\beta} \end{bmatrix}_{N+M} .
$$

Finally, the conditions (19a) to (19d) related to the estima-
tor Y^PCK(x_0) can thus be rewritten as follows, if the esti-
mation is made with a moving neighbourhood in which the
discontinuities vanishes thanks to the noises $n_\alpha$ and $n_\beta$:

$$ 
\begin{align*}
(K_{11}^{bb} + \mathbf{n}_{1b})\lambda_b + K_{12}^{bb}\nu_b + \mu_0 + \mu_j g_b^j &= K_{10}^b - (K_{11}^{bn} + \mathbf{n}_{1n})\bar{\lambda}_n - K_{12}^{bn}\bar{\nu}_n \\
(K_{11}^{nb} + \mathbf{n}_{1n})\lambda_b + K_{12}^{nb}\nu_b + \mu_0 + \mu_j g_n^j + \mathbf{1}_{11}^{nn}\tau_n &= K_{10}^n - (K_{11}^{nn}\bar{\lambda}_n + \mathbf{n}_{1n}) - K_{12}^{nn}\bar{\nu}_n \\
K_{21}^{bb}\lambda_b + (K_{22}^{bb}\nu_b + \mathbf{n}_{2b}) + \nu_0 + \nu_i f_b^j &= K_{20}^b - K_{21}^{bn}\bar{\lambda}_n - (K_{22}^{bn} + \mathbf{n}_{2n})\bar{\nu}_n \\
K_{21}^{nb}\lambda_b + (K_{22}^{nb} + \mathbf{n}_{2n})\nu_b + \nu_0 + \nu_i f_n^j + \mathbf{1}_{22}^{nn}\omega_n &= K_{20}^n - K_{21}^{nn}\bar{\lambda}_n - (K_{22}^{nn} + \mathbf{n}_{2n})\bar{\nu}_n
\end{align*}
$$

$\mathbf{n}_{1b}$ is the vector of noises $\{n(\mathbf{x}_\alpha)\}, \alpha \in \mathcal{B}_\alpha$.

$\mathbf{n}_{1n}$ is the vector of noises $\{n(\mathbf{x}_\alpha)\}, \alpha \in \bar{\mathcal{B}}_\alpha$.

$\mathbf{n}_{2b}$ is the vector of noises $\{n(\mathbf{x}_\beta)\}, \beta \in \mathcal{B}_\beta$.

$\mathbf{n}_{2n}$ is the vector of noises $\{n(\mathbf{x}_\beta)\}, \beta \in \bar{\mathcal{B}}_\beta$.

Let us note that in air quality kriging-based maps, the discontinuities mainly arises on very large domain of estimations, over Europe for instance, where the monitoring network is not homogeneously distributed: a lot of data are available in Western and Central Europe while the network is sparse elsewhere. Thus, to enable a local fitting of the drift, a number-based neighborhood $\delta_{\alpha 0} = \{\alpha, ||\mathbf{x}_\alpha - \mathbf{x}_0|| < ||\mathbf{x}_N - \mathbf{x}_0||\}$, where $\mathbf{x}_N$ denotes the Nth nearest neighbour of $\mathbf{x}_0$, is preferred to a distance-based neighborhood $\delta_{\alpha 0} = \{\alpha, ||\mathbf{x}_\alpha - \mathbf{x}_0|| \le D\}$; the latter including too many stations where the network is dense, and too few in the badly informed areas. In Rivoirard and Romary [31], the continuous kriging is distance-based driven. To overcome this problem, and because the estimations are done on regular grids, the distance $R(\mathbf{x}_0)$ between the target location $\mathbf{x}_0$ and the Nth nearest neighbor of $\mathbf{x}_0$ is computed for each gridcell. This distance spatially varies but in a continuous way since the grid is regular (see Fig. B.6). As a consequence, the distance-based neighborhood approach is kept but the radius used is no longer spatially constant: $R(\mathbf{x}_0)$ is substituted to $R$ in Eq.(26).

<figure><img src="image_8.png" /><figcaption>Figure B.6: Distance between the gricells $\mathbf{x}_0$ and the 20th nearest neighbour</figcaption></figure>
**Declaration of interests**

[x] The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[ ] The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

<figure><img src="image_9.png" /></figure>