Optimizing Air Pollution Modeling with a Highly-Convergent Quasi-Monte Carlo Method: A Case Study on the UNI-DEM Framework

Abstract

In this study, we present the development of an advanced air pollution modeling approach, which incorporates cutting-edge stochastic techniques for large-scale simulations of long-range air pollutant transportation. The Unified Danish Eulerian Model (UNI-DEM) serves as a crucial mathematical framework with numerous applications in studies concerning the detrimental effects of heightened air pollution levels. We employ the UNI-DEM model in our research to obtain trustworthy insights into critical questions pertaining to environmental preservation. Our proposed methodology is a highly convergent quasi-Monte Carlo technique that relies on a unique symmetrization lattice rule. By fusing the concepts of special functions and optimal generating vectors, we create a novel algorithm grounded in the component-by-component construction method, which has been recently introduced. This amalgamation yields particularly impressive outcomes for lower-dimensional cases, substantially enhancing the performance of the most advanced existing methods for calculating the Sobol sensitivity indices of the UNI-DEM model. This improvement is vital, as these indices form an essential component of the digital ecosystem for environmental analysis.

1. Introduction

Environmental security is a major global concern and poses many challenges due to its sensitive nature for society and healthcare systems. To ensure the accuracy and reliability of large-scale air pollution computational models, multidimensional sensitivity analysis plays a crucial role in the validation process. Sensitivity analysis (SA) is the tool for checking how uncertainty in input data affects the accuracy of output results [1,2,3,4,5]. When modelling a complex phenomenon, multidimensional sensitivity analysis (MSA) can be a challenging task [6,7,8,9]. The main issue in sensitivity analysis is evaluating total sensitivity indices, which involves calculating a set of multiple integrals [10,11,12,13]. Monte Carlo methods are the most effective approach for solving these integrals [10,14,15,16,17,18].

The input data for sensitivity analysis is derived from the Unified Danish Eulerian Model (UNI-DEM) [19,20,21,22,23], which covers a large geographical region and is a fundamental tool for creating important digital twins, called Digital Air, developed recently in [24]. It is shown that the digital ecosystem consists of many numerical algorithms, some splitting techniques [25], many graphical tools, some useful scenarios, lots of meteorological and emission data [26], and a huge amount of geographical information. In this paper, we aim to improve the numerical algorithms for the digital twin.

The model can be defined using a function that represents it, as per its definition [27]:

$$\mathrm{u}=f\left(\mathrm{x}\right),\mathrm{x}=(x_1,x_2,\dots ,x_d)\in U^d\equiv {[0;1]}^d.$$

The Sobol approach is based on representing $f\left(\mathrm{x}\right)$ with a constant $f_0$ in the following manner [27]:

$$f\left(\mathrm{x}\right)=f_0+\sum _{\nu =1}^d\sum _{l_1<\dots <l_{\nu }}{f}_{l_1\dots l_{\nu }}(x_{l_1},x_{l_2},\dots ,x_{l_{\nu }}).$$

The expression (1) is commonly referred to as the ANOVA representation of $f\left(\mathrm{x}\right)$ according to [27]:

$${\int }_0^1{f}_{l_1\dots l_{\nu }}(x_{l_1},x_{l_2},\dots ,x_{l_{\nu }})\mathrm{d}{x}_{l_k}=0,1\le k\le \nu ,\nu =1,\cdots ,d.$$

The metrics

$${\displaystyle \mathbf{D}={\int }_{U^d}{f}^2\left(\mathrm{x}\right)\mathrm{d}\mathrm{x}-f_0^2,{\mathbf{D}}_{l_1\dots l_{\nu }}=\int f_{l_1\dots l_{\nu }}^2\mathrm{d}{x}_{l_1}\dots \mathrm{d}{x}_{l_{\nu }}}$$

are named total and partial variances [27]. The total variance is as follows:

$$\mathbf{D}=\sum _{\nu =1}^d\sum _{l_1<\dots <l_{\nu }}{\mathbf{D}}_{l_1\dots l_{\nu }}.$$

The Sobol global sensitivity indices (GSIs) [4,27] are defined by

$${\displaystyle S_{l_1\dots l_{\nu }}=\frac{{\mathbf{D}}_{l_1\dots l_{\nu }}}{\mathbf{D}},\nu \in \{1,\dots ,d\}.}$$

Thus, the total sensitivity index (TSI) of input parameter $x_i,i\in \{1,\dots ,d\}$ is calculated according to [28]:

$${\displaystyle TSI\left(x_i\right)=S_i+\sum _{l_1\ne i}{S}_{i{l}_1}+\sum _{l_1,l_2\ne i,l_1<l_2}{S}_{i{l}_1{l}_2}+\dots +S_{i{l}_1\dots l_{d-1}},}$$

where $S_{i{l}_1\dots l_{j-1}}$ is the j^-th order SI for $x_i(2\le j\le d)$.

As per [29], the term $S_i$ is referred to as the “main effect” of $x_i$ when $j=1$. In the case of $j=2$, $S_{ij}$ is denoted as “two-way interactions” (second-order SIs). Similarly, for $j=3$, $S_{ijk}$ is termed as “three-way interactions” (third-order SIs), etc. In this study we will be interested in all important GSIs of first and second order, especially those that are small in value.

The total sensitivity of the output variance to an input parameter $x_i,i\in \{1,\dots ,d\}$ is represented as [30]:

$${\displaystyle S_{T_i}=S_i+\sum _{j\ne i}{S}_{ij}+\sum _{j\ne i,k\ne i,j<k}{S}_{ijk}+\dots .}$$

This demonstrates that applying Sobol’s approach for MSA transforms it into a task involving multidimensional numerical integration [31].

2. Methods and Algorithms

Let us examine the given objective for performing integration in multiple dimensions, specifically in dimension s:

$$I\left(f\right):=I={\int }_{U^s}f\left(x\right)\mathrm{d}\mathrm{x}.$$

(1)

Lattice rules utilize deterministic sequences instead of random sequences, distinguishing them as a unique form of low-discrepancy sequences. These sequences possess a specific structure that ensures a more even distribution across the domain. It has been observed that when the integral exhibits sufficient regularity, lattice rules tend to surpass not only conventional Monte Carlo methods but also various other low-discrepancy sequences in terms of performance and accuracy. The inherent advantages of lattice rules stem from their ability to efficiently capture the underlying patterns and optimize the integration process, making them a valuable tool in numerical computations.

Let us consider the quadrature formula

$$I_N\left(f\right)=\frac{1}{N}\sum _{i=1}^{N}f\left(x_i\right),$$

(2)

where $P_N=x_1,x_2,\cdots ,x_N$ represents the integration nodes for (2), where $x_i\in {[0,1)}^s$. We will utilize the integration nodes from (2) as provided in [32,33]:

$$x_k=\left(\left\{\frac{k{z}_1}{N}\right\},\left\{\frac{k{z}_2}{N}\right\},\cdots ,\left\{\frac{k{z}_s}{N}\right\}\right),k=1,2,\cdots ,N.$$

(3)

Here, N represents the count of nodes or points, z denotes a generating vector with dimensions s for the lattice set, and a refers to the fractional component of a, obtained by $\left\{a\right\}=a-\left[a\right]$. Now Formula (2) with nodes (3) and generators z are called rank-1 lattice rules [34]. We will use a special type of rank-1 lattice (R1L).

Now we propose a special R1L, defined in the following way. In the one-dimensional case, we define a function appropriate for integrand functions that exhibit periodicity to be utilized in conjunction with a non-periodic function F by applying the R1L to the function

$$L\left(x\right)=\left(F\left(x\right)+F(1-x)\right)/2,$$

in one dimension. Considering two dimensions, the function L is established as

$$L(x_1,x_2)=\left(F(x_1,x_2)+F(x_1,1-x_2)+F(1-x_1,x_2)+F(1-x_1,1-x_2)\right)/4.$$

The function $L(x_1,\dots ,x_s)$ is extended to accommodate s dimensions in the given definition:

$$L(x_1,\dots ,x_s)=2^{-s}\sum _{\epsilon \in {\{0,1\}}^s}F\left({\epsilon }_1{x}_1+(1-{\epsilon }_1)(1-x_1),\dots ,{\epsilon }_s{x}_s+(1-{\epsilon }_s)(1-x_s)\right).$$

(4)

The terms being summed can be viewed as vertices of a parallelotope, with the diagonals of the parallelotope intersecting at a specific point: $\left(1/2,1/2,\dots ,1/2\right)\in {[0,1]}^s$. Formula (4) is equivalent to

$$L(x_1,\dots ,x_s)=\sum _{\epsilon \in {\{0,1\}}^s}F\left(x_1^{{\epsilon }_1}{(1-x_1)}^{1-{\epsilon }_1},x_2^{{\epsilon }_2}{(1-x_2)}^{1-{\epsilon }_2},\dots ,x_s^{{\epsilon }_s}{(1-x_s)}^{1-{\epsilon }_s}\right).$$

Finding the most suitable generating vectors in high-dimensional spaces with a fixed value of N presents significant computational challenges. Selecting an effective generator vector is the primary obstacle in R1L. It is preferable for the one-dimensional projections of R1L to have N unique values, necessitating restrictions on each component of $\mathbf{z}$ to satisfy this requirement:

$$U_N:=\left\{z\in \mathbb{N}:1\le z\le N-1,\gcd (z,N)=1\right\}.$$

Actually, the size of $U_N$ is $\left|U_N\right|=\phi \left(N\right)$, which is the Euler totient function. If the prime factorization of N is defined as $N=p_1^{{\alpha }_1}{p}_2^{{\alpha }_2}\cdots p_k^{{\alpha }_k}$, it is true that

$$\phi \left(N\right)=(p_1^{{\alpha }_1}-p_1^{{\alpha }_1-1})(p_2^{{\alpha }_2}-p_2^{{\alpha }_2-1})\cdots (p_k^{{\alpha }_k}-p_k^{{\alpha }_k-1}).$$

As N becomes large, the function $\phi \left(N\right)$ increases approximately at a rate proportional to N, specifically expressed as $\frac{1}{\phi \left(N\right)}=\mathcal{O}\left(\frac{loglogN}{N}\right)$. When N is a prime number, $\phi \left(N\right)$ is equal to $N-1$. Consequently, there are $N-1$ potential options for each component of $\mathbf{z}$, and a total of ${(N-1)}^s$ possibilities for the generating vector $\mathbf{z}$.

If N and s are both large, it is practically impossible to use brute force to search for the generating vector that meets a predetermined error criterion. Therefore, advanced number theory methods have been developed, including those based on criteria such as the Zaremba index and worst function errors. One example of a plain construction algorithm is the Korobov algorithm [33].

If an integer a is chosen in a manner that $1\le a\le N-1$ and $\gcd (a,N)=1$, then

$$\mathbf{z}=\mathbf{z}\left(a\right):=(1,a,a^2,a^3,\dots ,a^{s-1})\mathrm{mod}N.$$

It is clear that the parameter a has up to $N-1$ options, resulting in the same number of choices for the generating vector $\mathbf{z}$. In some cases, it may be feasible to check all $N-1$ choices for $\mathbf{z}$ and choose the one with the best properties. However, when N and s are relatively large, this brute-force approach becomes impractical. Instead, we can use a more sophisticated technique called the component-by-component (CBC) construction method, which has been recently developed and described by Nuyens and Kuo [35]. This method involves obtaining different generating vectors through a component-wise construction process.

Initially, the value of $z_1$ is set to 1. Subsequently, $z_1$ remains fixed, and a value for $z_2$ is selected from the set $U^N$, which consists of natural numbers z satisfying $1\le z\le N-1$ and $\gcd (z,N)=1$. The choice of $z_2$ aims to minimize the predefined error criterion in the two-dimensional case, as described in [35]. This process is then repeated iteratively for $i=3,\dots ,s$, with $z_i$ chosen from $U^N$ in a manner that minimizes the predefined error criterion in i dimensions.

We will now generate a set of special generators (SGs) for the R1L method. The first SG, referred to as 1PT, will be created using the CBC construction of a rank-1 lattice rule with a prime number of points and product weights. The second SG, called 1OD, will also use the CBC construction of a rank-1 lattice rule with a prime number of points, but with order-dependent weights. For the third SG, denoted as 1EXPT, the CBC construction of a rank-1 lattice sequence with a prime power of points and product weights will be employed. Lastly, the fourth SG, 1EXOD, will utilize the CBC construction of a rank-1 lattice sequence with a prime power of points and order-dependent weights.

The fifth SG will be the special polynomial construction following the paradigm of rank-1 LS in base 2 with product weights, denoted by 1POLY. On the first stage of R1L,

$$\mathbf{z}=(z_1,z_2,\dots z_s)$$

is obtained by the CBC. So, the lattice points are built following

$${\mathbf{x}}_k=\left\{\frac{k}{N}\mathbf{z}\right\},k=1,\cdots ,N.$$

(5)

Finally, the approximate value $I_N$ of the integral is derived:

$$I_N=\frac{1}{N}\sum _{k=1}^{N}f\left(\left\{\frac{k}{N}\mathbf{z}\right\}\right).$$

(6)

Lattice sets prove to be highly suitable for integrands characterized by their smoothness. The key benefit of employing the lattice method lies in its linear computational complexity, which leads to significant time savings when evaluating multidimensional integrals. Notably, the number of calculations needed to generate the generating vector is asymptotically lower than $\mathcal{O}\left(N\right)$. Furthermore, the generation of each new point within the lattice set necessitates a constant number of operations. As a result, to obtain a lattice set consisting of N points following the described procedure, only $\mathcal{O}\left(N\right)$ operations are required.

A method is shown for constructing effective lattice rules, utilizing the fast component-by-component algorithm. This algorithm enables the construction of high-quality lattice rules with a time complexity of $\mathcal{O}(sNlog(N\left)\right)$ and a memory requirement of $\mathcal{O}\left(N\right)$. The algorithm’s effectiveness could be demonstrated in the case where N is a prime number and the underlying function space is a weighted, shift-invariant, and tensor-product reproducing kernel Hilbert space [36,37]. Remarkably, it is established that with a minor increase in construction cost, the fast algorithm can also be applied to more general weighted function spaces. Specifically, it proves advantageous for handling order-dependent weights. This extension broadens the applicability of the fast algorithm, enhancing its utility in various scenarios. In the next section, if $n=2^p$, we use the lowest prime $N\ge n$.

To benchmark the aforementioned approaches, we compare the performance of the lattice sequences with the most basic Monte Carlo approach—the crude Monte Carlo—as well as test against two of the good quasi-Monte Carlo methods, the Halton and Sobol sequences. They are both types of quasi-Monte Carlo sequences used in numerical integration and simulation. While they share similarities in their underlying principles, there are notable differences between them. To begin with, the Halton sequence is constructed using a specific base for each dimension. Each dimension in the sequence corresponds to a prime number as its base. The sequence is created by generating a series of points in the unit interval, with each point having a unique combination of prime base fractions. The Sobol sequence, however, is constructed based on the concept of Gray codes and primitive polynomials. It utilizes a set of carefully chosen irreducible polynomials to generate unique points in each dimension. The Sobol sequence ensures a more even distribution of points in the unit interval compared to traditional random sequences.

What is more, the Halton sequence exhibits relatively good discrepancy properties, especially for low-dimensional problems. However, as the number of dimensions increases, the Halton sequence may suffer from more pronounced clustering. The Sobol sequence is designed to have superior discrepancy properties, even for high-dimensional problems. It distributes points more uniformly across the integration space, reducing the clustering effect observed in higher dimensions.

Moreover, the Halton sequence is relatively straightforward to implement and requires minimal computational resources. It is commonly used for low-dimensional problems and situations where simplicity and efficiency are valued. In contrast, the Sobol sequence can be more complex to implement due to the need for primitive polynomials. However, it offers greater flexibility and superior performance for high-dimensional problems, making it suitable for more demanding applications. Furthermore, the Sobol sequences allow for different optimizations, while the Halton sequences cannot be seriously modified.

3. Results and Discussion

In this study, the advanced stochastic lattice algorithms described in the previous section, R1L-1PT, R1L-1OD, R1L-1EXPT, R1L-1EXOD, and R1L-1POLY, are compared with the best available method for estimating Sobol indices, namely the Sobol sequence (Sobol) [38,39,40] and Halton sequence (Halton) [41,42], and the basic Monte Carlo algorithm (Crude) [18]. We employ two approaches to conduct sensitivity studies regarding emission levels (referred to as SSREL) and sensitivity studies concerning specific chemical reaction rates (referred to as SSRCRR) in order to analyze the concentration variations in pollutants in UNI-DEM [43,44]. The quantity being estimated is denoted as EQ, with its corresponding reference value indicated as RV. The relative error is represented by RE, and an approximate evaluation is denoted by AE.

Chemical reactions have a crucial significance within the UNI-DEM framework. The mathematical model relies on the interconnection of equations through these chemical reactions. Furthermore, the system of equations exhibits nonlinearity and stiffness primarily due to the involvement of chemical processes. The specific chemical scheme employed in the model is the condensed CBM-IV (Carbon Bond Mechanism). In its current state, additional reactions have been incorporated to account for the formation of intricate positive ammonium ions in the atmosphere. This involves nitrogen atoms from ammonia molecules bonding with hydrogen ions generated through the dissociation of water, utilizing their available electron pairs.

The intensity of each reaction is dictated by its corresponding chemical rate constant, which can be regarded as a stochastic variable for the purpose of sensitivity analysis. In the initial phase of our investigation, the objective is to identify the key chemical rate constants based on a specific criterion. Ozone, a highly hazardous air pollutant, is chosen as the primary criterion, specifically focusing on mean monthly ozone concentration. Through extensive experimentation involving the perturbation of numerous coefficients associated with suspected chemical reactions, it is revealed f. i. that chemical reaction #22, $H{O}_2+NO\to OH+N{O}_2$ (time-dependent), is one of the utmost critical reactions.

Considering SSREL, we will investigate the sensitivity of the model output with regards to the variation of input emissions of anthropogenic pollutants. Specifically, we will be looking at the mean monthly concentrations of various critical pollutants, with ammonia in Milan being a particular area of interest. Our goal is to determine how sensitive the model is to changes in these input emissions, which can help us better understand the accuracy and reliability of our predictions.

The input emissions of anthropogenic pollutants consists of four distinct components, and we will analyze their effects on the model output. The components $\mathbf{E}=({\mathbf{E}}^{\mathbf{A}},{\mathbf{E}}^{\mathbf{N}},{\mathbf{E}}^{\mathbf{S}},{\mathbf{E}}^{\mathbf{C}})$ themselves are defined as follows:

$$\begin{array}{cc}{\mathbf{E}}^{\mathbf{A}}-\mathrm{ammonia}\left(N{H}_3\right);\hfill & {\mathbf{E}}^{\mathbf{S}}-\mathrm{sulphur} \mathrm{dioxide}\left(S{O}_2\right);\hfill \\ {\mathbf{E}}^{\mathbf{N}}-\mathrm{nitrogen} \mathrm{oxides}(NO+N{O}_2);\hfill & {\mathbf{E}}^{\mathbf{C}}-\mathrm{anthropogenic} \mathrm{hydrocarbons}.\hfill \end{array}$$

The output of the model is the average monthly concentration of the following three pollutants:

• $s_1$—ozone ($O_3$);
• $s_2$—ammonia ($N{H}_3$);
• $s_3$—ammonium sulfate and ammonium nitrate ($N{H}_{4}S{O}_4+N{H}_{4}N{O}_3$).

The results for REs for the AE of the $f_0$, $\mathbf{D}$, $S_i$, and $S_i^{\mathrm{tot}}$ are obtained using Crude, Sobol, Halton, R1L-1PT, R1L-1OD, R1L-1EXPT, R1L-1EXOD, R1L-1POLY. The results for different quantities are presented in

Table 1. Relative error of approximating the model function $f_0\approx 0.048$.

n	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$2^4$	7.9e-02	2.0e-02	1.0e-03	6.4e-04	1.1e-04	9.5e-05	9.5e-05	5.9e-04
$2^6$	1.7e-02	5.0e-03	2.5e-03	1.4e-04	7.1e-06	5.9e-06	5.9e-06	3.7e-05
$2^8$	1.7e-02	1.2e-03	6.7e-04	8.9e-06	4.5e-07	3.7e-07	3.7e-07	2.3e-06
$2^{10}$	1.2e-03	3.1e-04	2.5e-04	5.6e-07	2.8e-08	2.3e-08	2.3e-08	1.4e-07
$2^{12}$	3.1e-03	7.8e-05	1.4e-04	3.5e-08	1.7e-09	1.5e-09	1.5e-09	9.0e-09
$2^{14}$	3.8e-03	2.0e-05	2.8e-05	2.2e-09	1.1e-10	9.1e-11	9.1e-11	5.6e-10

,

Table 2. Relative error of approximating the total variance $\mathbf{D}\approx 0.0002$.

n	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$2^4$	1.8e+00	6.1e-02	1.3e+00	4.4e-02	7.2e-03	7.9e-03	6.9e-03	2.2e-01
$2^6$	5.1e-01	1.2e-02	2.2e-01	1.7e-03	4.7e-04	4.7e-04	4.9e-04	1.4e-02
$2^8$	2.1e-01	2.5e-03	9.3e-02	1.2e-04	4.8e-07	3.1e-05	3.1e-05	8.5e-04
$2^{10}$	4.4e-02	4.8e-05	1.5e-02	5.5e-06	1.2e-08	1.7e-06	1.5e-06	5.3e-05
$2^{12}$	1.3e-01	3.4e-05	8.1e-03	4.3e-07	1.1e-08	1.2e-07	1.1e-07	3.4e-06
$2^{14}$	9.9e-02	1.1e-05	1.9e-03	4.5e-08	3.8e-08	6.0e-09	1.0e-08	3.4e-07

,

Table 3. Relative error of approximating the sensitivity indices and the total sensitivity indices for $n=2^6$.

SI	EQ	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
		Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$S_1$	9e-01	6.9e-01	5.8e-03	8.5e-02	1.9e-03	1.8e-04	2.3e-05	1.4e-05	1.3e-02
$S_2$	2e-04	5.9e+01	1.4e+01	6.6e+00	6.3e-02	2.7e-01	2.3e-03	6.2e-03	4.9e-02
$S_3$	1e-01	6.1e+00	1.4e-02	7.7e-01	2.2e-03	1.7e-03	1.5e-05	8.2e-06	1.1e-01
$S_4$	4e-05	2.7e+02	9.8e-01	2.4e+01	3.6e+01	1.1e-01	2.1e-02	1.5e-02	4.8e-02
$S_1^{\mathrm{tot}}$	9e-01	7.5e-01	2.2e-04	9.8e-02	1.9e-03	1.5e-04	1.9e-06	4.6e-07	1.3e-02
$S_2^{\mathrm{tot}}$	2e-04	1.2e+02	9.0e-01	2.3e+00	4.0e-02	2.3e-01	5.1e-03	1.4e-03	4.6e-02
$S_3^{\mathrm{tot}}$	1e-01	5.7e+00	4.8e-02	6.6e-01	2.5e-03	1.9e-03	1.9e-04	1.2e-04	1.1e-01
$S_4^{\mathrm{tot}}$	5e-05	2.7e+02	2.6e+01	3.5e+01	3.0e+01	9.6e-02	1.2e-02	1.0e-03	3.7e-02

,

Table 4. Relative error of approximating the sensitivity indices and the total sensitivity indices for $n=2^{10}$.

SI	EQ	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
		Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$S_1$	9e-01	6.2e-02	3.6e-05	1.7e-02	8.8e-06	1.1e-06	6.9e-08	1.1e-07	5.2e-05
$S_2$	2e-04	4.4e+00	7.9e-02	4.0e-02	1.9e-04	1.0e-03	4.2e-05	1.6e-03	4.4e-05
$S_3$	1e-01	4.9e-01	3.3e-03	1.4e-01	2.0e-05	8.4e-06	1.8e-07	3.7e-07	4.2e-04
$S_4$	4e-05	5.4e+00	1.7e+00	1.9e+00	1.4e-01	6.3e-04	3.6e-03	4.8e-05	1.6e-04
$S_1^{\mathrm{tot}}$	9e-01	6.0e-02	3.6e-04	1.7e-02	8.7e-06	8.3e-07	1.5e-07	3.4e-07	5.2e-05
$S_2^{\mathrm{tot}}$	2e-04	1.7e-01	1.3e-01	7.4e-01	1.0e-04	8.5e-04	1.8e-05	1.4e-03	1.1e-04
$S_3^{\mathrm{tot}}$	1e-01	5.0e-01	4.8e-04	1.4e-01	2.0e-05	1.0e-05	1.9e-06	2.2e-06	4.2e-04
$S_4^{\mathrm{tot}}$	5e-05	2.9e+00	7.2e-04	2.4e+00	1.2e-01	5.9e-04	2.9e-03	1.8e-05	3.3e-04

,

Table 5. Relative error of approximating the sensitivity indices and the total sensitivity indices for $n=2^{14}$.

SI	EQ	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
		Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$S_1$	9e-01	2.4e-02	2.6e-06	1.2e-03	3.1e-08	1.3e-09	9.7e-10	4.0e-09	2.7e-07
$S_2$	2e-04	1.3e-01	1.4e-03	3.4e-02	2.6e-06	5.4e-06	8.8e-07	3.1e-07	3.9e-07
$S_3$	1e-01	1.8e-01	1.1e-05	1.0e-02	5.3e-08	5.2e-09	1.7e-09	1.8e-09	2.2e-06
$S_4$	4e-05	7.4e+00	9.6e-03	4.0e-01	5.4e-04	8.9e-06	1.2e-06	3.0e-06	1.3e-06
$S_1^{\mathrm{tot}}$	9e-01	2.1e-02	1.9e-06	1.3e-03	3.0e-08	2.2e-09	3.8e-11	1.5e-10	2.7e-07
$S_2^{\mathrm{tot}}$	2e-04	3.3e+00	6.6e-03	6.5e-02	1.1e-06	4.3e-06	6.7e-07	9.2e-07	4.8e-07
$S_3^{\mathrm{tot}}$	1e-01	2.0e-01	1.2e-05	9.7e-03	5.7e-08	2.5e-09	9.7e-09	2.9e-08	2.2e-06
$S_4^{\mathrm{tot}}$	5e-05	7.6e+00	2.5e-03	5.2e-01	4.5e-04	6.1e-06	8.0e-07	8.2e-06	9.6e-07

,

Table 6. Relative error of approximating the model function $f_0\approx 0.27$.

n	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$2^4$	4.4e-02	5.5e-03	8.3e-04	1.4e-02	1.0e-03	1.1e-04	1.1e-04	1.3e-02
$2^6$	1.1e-02	1.6e-03	2.0e-03	3.5e-03	6.3e-05	7.0e-06	7.0e-06	8.3e-04
$2^8$	1.2e-02	2.9e-04	6.3e-04	2.2e-04	3.9e-06	4.4e-07	4.4e-07	5.2e-05
$2^{10}$	6.8e-03	7.5e-05	3.7e-04	1.4e-05	2.5e-07	2.7e-08	2.7e-08	3.2e-06
$2^{12}$	7.2e-03	1.9e-05	8.6e-05	8.5e-07	1.5e-08	1.7e-09	1.7e-09	2.0e-07
$2^{14}$	4.2e-04	5.1e-06	3.4e-05	5.3e-08	9.6e-10	1.1e-10	1.1e-10	1.3e-08

,

Table 7. Relative error of approximating the total variance $\mathbf{D}\approx 0.0025$.

n	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$2^4$	1.6e+00	3.7e-01	9.3e-01	2.0e-01	1.4e-02	1.8e-02	8.3e-03	4.8e-01
$2^6$	1.0e+00	1.0e-02	3.3e-01	3.2e-01	4.9e-05	3.3e-04	3.5e-03	3.0e-02
$2^8$	3.3e-01	6.7e-03	2.7e-02	2.0e-02	1.8e-04	1.6e-04	3.8e-04	1.9e-03
$2^{10}$	3.0e-01	1.2e-03	5.0e-02	1.0e-03	3.8e-08	2.4e-05	1.3e-04	1.4e-04
$2^{12}$	9.2e-02	1.8e-04	2.8e-03	6.4e-05	2.0e-07	3.0e-06	8.2e-06	7.0e-06
$2^{14}$	5.9e-02	8.0e-05	2.4e-03	5.3e-06	4.7e-08	9.2e-08	2.9e-07	1.4e-06

,

Table 8. Relative error of approximating the sensitivity indices and the total sensitivity indices for $n=2^6$.

SI	EQ	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
		Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$S_1$	4e-01	9.6e-01	4.4e-02	8.9e-01	4.8e-01	3.8e-03	2.0e-04	3.1e-03	2.7e-02
$S_2$	3e-01	4.7e-02	4.8e-02	9.3e-01	4.7e-01	4.1e-03	5.7e-04	2.6e-03	2.9e-02
$S_3$	5e-02	2.1e+00	1.4e+00	1.7e+00	1.1e+01	1.1e-01	2.7e-04	1.6e-02	5.6e-01
$S_4$	3e-01	1.4e+00	2.9e-02	7.4e-01	6.8e-01	5.3e-03	8.5e-05	3.3e-03	2.9e-02
$S_5$	4e-07	9.2e+03	1.0e+03	7.6e+03	2.7e+01	1.9e+00	5.3e-02	1.0e-02	4.5e+00
$S_6$	2e-02	4.4e+00	2.0e+00	2.4e+00	5.4e-01	4.1e-02	5.1e-03	3.6e-04	3.0e-02
$S_1^{\mathrm{tot}}$	4e-01	1.3e+00	7.4e-02	6.8e-01	4.9e-01	4.2e-03	2.9e-04	2.6e-03	2.6e-02
$S_2^{\mathrm{tot}}$	3e-01	1.5e-02	3.6e-02	1.2e+00	4.8e-01	5.0e-03	2.7e-04	3.8e-03	2.8e-02
$S_3^{\mathrm{tot}}$	5e-02	1.4e+00	1.2e+00	1.4e+00	1.1e+01	1.1e-01	8.4e-04	1.6e-02	5.5e-01
$S_4^{\mathrm{tot}}$	3e-01	1.6e+00	6.9e-02	4.0e-01	6.9e-01	6.3e-03	7.2e-04	1.8e-03	2.8e-02
$S_5^{\mathrm{tot}}$	2e-04	3.6e+01	1.3e+00	1.4e+01	5.8e-01	1.5e-01	3.0e-02	3.0e-02	7.3e-02
$S_6^{\mathrm{tot}}$	2e-02	3.3e+00	1.6e+00	2.4e+00	5.5e-01	4.1e-02	4.3e-03	1.3e-03	2.8e-02
$S_{12}$	6e-03	1.9e+01	1.6e-01	5.3e+00	5.1e-01	1.3e-02	1.7e-02	1.4e-02	1.3e-02
$S_{14}$	5e-03	1.3e+00	4.3e+00	6.9e+00	1.2e+00	1.6e-02	1.4e-02	8.9e-03	1.7e-02
$S_{15}$	8e-06	2.3e+02	1.0e+01	3.7e+01	5.5e-01	2.9e-02	3.0e-02	1.1e-02	1.6e-02
$S_{24}$	3e-03	2.4e+01	5.3e+00	1.7e+01	1.2e+00	5.8e-02	4.0e-02	9.1e-02	4.2e-04
$S_{45}$	1e-05	4.8e+01	2.7e+01	5.5e+02	1.5e+00	4.7e-02	3.0e-02	8.9e-03	1.7e-01

,

Table 9. Relative error of approximating the sensitivity indices and the total sensitivity indices for $n=2^{10}$.

SI	EQ	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
		Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$S_1$	4e-01	1.6e-01	4.7e-03	8.9e-02	9.5e-04	2.3e-05	5.3e-06	1.2e-04	1.3e-04
$S_2$	3e-01	7.0e-02	1.7e-03	7.6e-02	1.0e-03	3.0e-05	1.8e-05	9.3e-05	1.0e-04
$S_3$	5e-02	3.1e-01	9.5e-03	2.2e-01	2.4e-02	4.7e-04	3.9e-05	1.6e-04	2.3e-03
$S_4$	3e-01	3.9e-01	1.1e-03	6.0e-02	1.6e-03	3.3e-05	2.4e-05	1.3e-04	1.6e-04
$S_5$	4e-07	1.5e+03	6.6e+02	1.1e+02	7.6e-02	2.1e-02	1.6e-02	1.7e-02	2.2e-02
$S_6$	2e-02	2.4e-01	6.0e-03	6.7e-01	1.0e-03	3.3e-04	7.2e-05	2.8e-03	3.7e-04
$S_1^{\mathrm{tot}}$	4e-01	2.6e-01	1.4e-03	7.6e-02	9.8e-04	3.5e-05	2.5e-07	1.1e-04	1.2e-04
$S_2^{\mathrm{tot}}$	3e-01	1.3e-01	2.5e-03	7.8e-02	1.0e-03	3.9e-05	2.1e-05	9.8e-05	8.4e-05
$S_3^{\mathrm{tot}}$	5e-02	1.6e-01	6.7e-04	2.6e-01	2.3e-02	4.4e-04	2.2e-05	1.6e-04	2.2e-03
$S_4^{\mathrm{tot}}$	3e-01	4.1e-01	4.8e-04	7.9e-02	1.6e-03	1.5e-05	1.4e-05	1.2e-04	1.3e-04
$S_5^{\mathrm{tot}}$	2e-04	3.3e+01	1.1e+00	6.1e+00	2.2e-03	1.5e-03	6.2e-05	8.2e-04	5.2e-04
$S_6^{\mathrm{tot}}$	2e-02	1.1e+00	5.4e-03	7.6e-01	1.1e-03	3.5e-04	7.8e-05	2.8e-03	3.8e-04
$S_{12}$	6e-03	5.2e+00	7.4e-02	1.0e+00	1.2e-03	2.4e-04	6.3e-05	4.2e-05	2.5e-05
$S_{14}$	5e-03	8.8e-02	7.3e-02	4.6e-01	3.1e-03	6.0e-04	3.0e-04	2.3e-05	4.9e-04
$S_{15}$	8e-06	6.6e+01	1.6e+01	2.4e+01	1.3e-03	1.4e-03	1.9e-04	1.9e-04	5.5e-05
$S_{24}$	3e-03	3.3e+00	3.4e-02	1.1e+00	3.0e-03	3.2e-04	1.6e-04	2.5e-06	1.3e-03
$S_{45}$	1e-05	3.1e+01	4.6e+00	4.5e+01	6.5e-03	3.7e-03	2.9e-04	1.6e-03	5.1e-05

and

Table 10. Relative error of approximating the sensitivity indices and the total sensitivity indices for $n=2^{14}$.

SI	EQ	Crude	Sobol	Halton	R1L-1PT	R1L-1OD	R1L-1EXPT	R1L-1EXOD	R1L-1POLY
		Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.	Rel.
$S_1$	4e-01	5.8e-02	4.4e-05	9.5e-03	5.2e-06	2.1e-07	9.6e-08	2.8e-07	4.2e-06
$S_2$	3e-01	2.8e-01	1.1e-04	1.2e-02	5.0e-06	1.0e-07	1.8e-08	3.4e-07	1.7e-06
$S_3$	5e-02	5.7e-01	1.4e-04	2.3e-02	1.1e-04	2.2e-06	2.9e-07	1.7e-07	1.1e-05
$S_4$	3e-01	2.6e-01	3.6e-04	7.2e-03	6.1e-06	1.4e-07	3.6e-07	4.3e-07	1.5e-06
$S_5$	4e-07	1.1e+03	1.7e+01	3.1e+02	6.4e-05	4.1e-03	2.1e-05	3.5e-05	4.9e-03
$S_6$	2e-02	9.6e-02	3.0e-03	1.7e-02	1.3e-05	6.8e-07	1.2e-06	4.5e-06	1.8e-06
$S_1^{\mathrm{tot}}$	4e-01	5.0e-02	1.4e-04	7.9e-03	5.5e-06	2.7e-07	1.2e-07	3.1e-07	4.0e-06
$S_2^{\mathrm{tot}}$	3e-01	2.9e-01	1.4e-04	1.1e-02	5.3e-06	2.7e-08	4.8e-08	2.3e-07	1.7e-06
$S_3^{\mathrm{tot}}$	5e-02	5.7e-01	8.5e-05	2.3e-02	1.0e-04	1.6e-06	3.1e-07	2.2e-07	1.1e-05
$S_4^{\mathrm{tot}}$	3e-01	2.9e-01	3.4e-04	6.4e-03	6.5e-06	6.5e-08	3.0e-07	2.6e-07	1.6e-06
$S_5^{\mathrm{tot}}$	2e-04	4.1e+00	4.1e-04	6.0e-02	1.5e-05	1.0e-04	2.3e-06	1.5e-06	1.7e-04
$S_6^{\mathrm{tot}}$	2e-02	2.3e-01	1.6e-04	2.4e-02	8.6e-06	7.9e-07	1.3e-06	3.7e-06	1.9e-06
$S_{12}$	6e-03	8.6e-01	2.5e-03	4.7e-02	1.5e-05	7.1e-07	7.0e-07	9.0e-07	1.8e-06
$S_{14}$	5e-03	1.3e+00	1.1e-02	1.5e-02	1.3e-05	3.3e-06	9.3e-07	1.1e-06	2.5e-06
$S_{15}$	8e-06	2.0e+01	4.9e-02	1.7e+00	2.0e-05	5.7e-07	1.7e-06	9.5e-07	4.9e-06
$S_{24}$	3e-03	3.2e-01	4.9e-04	3.0e-02	1.3e-05	9.0e-06	1.4e-06	4.8e-05	2.5e-06
$S_{45}$	1e-05	3.3e+00	6.5e-03	8.6e-01	6.8e-05	1.5e-05	3.9e-06	1.0e-06	9.7e-06

. Specifically in the case of SSREL, the quantity $f_0$ is evaluated using a four-dimensional integral, while the remaining quantities are calculated using eight-dimensional integrals. Similarly, the same methods are applied to determine the relative errors (Rel) and approximate evaluations for the SSRCRR. In this case, $f_0$ is evaluated using a six-dimensional integral, while the other quantities involve twelve-dimensional integrals. The smallest errors are denoted with bold.

In case of the SSREL, one may observe the following. First, it is obvious that the constructed algorithm significantly improves the results by the best available Sobol sequence for calculating the indices. The standard Monte Carlo plain algorithm suffers from a loss of accuracy, and the Halton sequence for most of the cases gives only mediocre results. In Table 1 for the model function $f_0$, the best algorithm for all numbers of samples are R1L-1EXPT and R1L-1EXOD and one should observe that they produce almost the same results; this is clearly visible in Figure 1. For the total variance $\mathbf{D}$, the best algorithm for almost all numbers of samples is R1L-1OD; see the results in Table 2. However, for $n=2^{14}$, L1R-1OD is slightly outperformed by, again, R1L-1EXPT and R1L-1EXOD; see Figure 1. The behavior of the algorithm can also be seen in Figure 1. From Table 3, Table 4 and Table 5, one can conclude that for all first-order sensitivity indices and total sensitivity indices, the best algorithms are again R1L-1EXPT and R1L-1EXOD. This could be explained by the fact that the generators, based on the CBC of the lattice sequence, are generally better than the generator based on the CBC of the lattice rules. R1L-1POLY yields the smallest relative error in only one of the cases and produces results very close to the other cases. R1L-1PT is generally outperformed by the other R1L in almost all of the cases. Some recent work [28] states that having the smallest possible sensitivity indices is the most important part of UNI-DEM. The relative errors for the sensitive indices smallest in value can be seen in Figure 2. The performance of the algorithms in the case of the sensitivity indices smallest in value is the following. For $S_2$, the best algorithm is R1L-1EXOD, for $S_4$, the best is R1L-1EXPT, for $S_2^{\mathrm{tot}}$, the best is R1L-2POLY, and for $S_4^{\mathrm{tot}}$, the best is R1L-1EXPT.

In the case of the SSRCRR, one may observe the following. To begin with, it is obvious that the constructed algorithm significantly improves the results of the best available Sobol sequence for calculating the indices, but here the difference is not so pronounced as in the lower-dimensional case, since the Sobol sequence is generally better for higher dimensions. The standard Monte Carlo approach is again the least accurate method, and the Halton sequence for most of the cases gives poor results. In Table 6 for the model function $f_0$, the best algorithms for all numbers of samples are R1L-1EXPT and R1L-1EXOD, and one should observe that the results again are very similar to the behavior in the previous case; this is clearly visible in Figure 3. For the total variance D, the best algorithm for almost all numbers of samples is R1L-1OD, as in the case of SSREL; see the results in Table 7 and Figure 3. From Table 8, Table 9 and Table 10, it could be concluded that for all first-order sensitivity indices and total sensitivity indices, the best algorithms are again R1L-1EXPT and R1L-1EXOD. But the difference here is that for the largest number of samples in Table 10, R1L-1OPT outperformed the other two algorithms. It is an interesting property of the generator with order-dependent weights for higher dimensions and needs further analysis, which is not trivial. Also, the generators, based on the CBC of the lattice sequences, are generally better than the generators, based on the CBC of the lattice rules. R1L-1POLY yields the smallest relative error in only one of the cases and produces results very close to the other cases. R1L-1PT is again generally outperformed by the other R1L in almost all of the cases. The relative errors for the sensitive indices smallest in value can be seen in Figure 4. The performance of the algorithms in the case of the sensitivity indices smallest in value is the following. For $S_5$, the best algorithm is R1L-1EXPT, for $S_5^{\mathrm{tot}}$, the best is R1L-1EXOD, for $S_{15}$, the best is R1L-1OD, and for the $S_{45}$, the best is R1L-1EXOD.

The general conclusion is that the proposed five new algorithms produce especially good results for lower dimensions and significantly improve the results for the best available approaches, even for higher-dimensional cases. However, there are better variations of the Sobol sequence, based on shifted Sobol points and improved Sobol generators like Joe-Kuo and Broda, and we need to compare our methods with these improved approaches. Furthermore, other super convergent methods exist, based on digital sequences with Sobol matrices, and this also will be explored in some further investigations.

4. Conclusions

In this paper we develop a new highly convergent lattice symmetrization rule. The novel lattice stochastic method significantly improves the existing algorithms for a particular digital twin, known as digital air. The new algorithm is combined with the recently developed component-by-component construction technique. The newly proposed lattice algorithms give dramatic improvement for lower dimensions and significantly improve the results for the best available approaches, even for higher-dimensional cases. The proposed lattice algorithms show particular improvement for the sensitivity indices smallest in value of UNI-DEM, which is of key importance for improving the numerical algorithms for the digital ecosystem. We also have shown that for some considerations, sensitivity indices small in value are important. To summarize, we successfully improved the digital ecosystem by improving the existing stochastic approaches by a specially designed lattice symmetrization approach. In future work, a deeper comparison will be made between the developed lattice approach and the best available digital sequences and best available Sobol sequence generators. The enhancements we have made to compute the sensitivity indices of the model will have a significant impact on the precision of evaluating agricultural losses. This, in turn, will enable a more accurate estimation of the consequences of detrimental emissions on human health. The importance of this cannot be overstated, as it will provide valuable insights into the potential health risks posed by these emissions. By using our improved methodology, we can better understand the link between emissions and their effects on both agriculture and human health, thereby allowing us to take more informed actions to mitigate their impact.