## Abstract

This paper presents an adaptive Kriging based method to perform uncertainty quantification (UQ) of the photoelectron sheath and dust levitation on the lunar surface. The objective of this study is to identify the upper and lower bounds of the electric potential and that of dust levitation height, given the intervals of model parameters in the one-dimensional (1D) photoelectron sheath model. To improve the calculation efficiency, we employ the widely used adaptive Kriging method (AKM). A task-oriented learning function and a stopping criterion are developed to train the Kriging model and customize the AKM. Experiment analysis shows that the proposed AKM is both accurate and efficient.

## 1 Introduction

The Moon is directly exposed to solar radiation and solar wind plasma (drifting protons and electrons) lacking an atmosphere and a global magnetic field. Consequently, the lunar surface is electrically charged by the bombardment of solar wind plasma and emission/collection of photoelectrons. Near the illuminated lunar surface, the plasma sheath is dominated by photoelectrons, thus usually referred to as “photoelectron sheath.” Additionally, dust grains on the lunar surface may get charged and levitated from the surface under the influence of the electric field within the plasma sheath as well as gravity. This work is motivated by the high computational cost associated with uncertainty quantification (UQ) analysis of plasma simulations using high-fidelity kinetic models such as particle-in-cell. The main quantities of interest (QoI) of this study are the vertical structure of the photoelectron sheath and its effects on levitation of dust grains with different sizes and electric charges.

Both the electric potential ($\varphi $) and the electric field ($E$) on lunar surface are determined by many parameters, such as solar wind drifting velocity ($vd$), electron temperature ($Te$), photoelectron temperature ($Tp$), density of ions at infinity ($ni,\u221e$), and density of photoelectrons ($np$). Due to uncertain factors in lunar environment, the electric potential, electric field, and the dust levitation height, etc., are also uncertain. While many sources uncertainty may exist, they are generally categorized as either aleatory or epistemic. Uncertainties are characterized as epistemic, if the modeler sees a possibility to reduce them by gathering more data or by refining models. Uncertainties are categorized as aleatory if the modeler does not foresee the possibility of reducing them [1]. An example of the aleatory uncertainty in lunar environment is the solar wind parameters, and an example of the epistemic uncertainty is the photoelectron temperature which is obtained by limited measurement data from Apollo missions. For lunar landing missions, one needs to take into consideration the uncertainties of the electrostatic and dust environment near the lunar surface. For example, the upper and lower bounds of the electric field and dust grain levitation heights in the photoelectron sheath should be considered when determining whether it is safe for a certain area to land a spacecraft.

Determining the bounds of the electric potential, electric field, and dust levitation height, however, is computationally expensive, because the particle-based kinetic models such as particle-in-cell simulations are time-consuming to evaluate. To address this issue, we develop an adaptive Kriging method (AKM), which can determine those bounds with a small number of calculations of the model. It is straightforward to train and obtain an accurate Kriging model [2] to replace the actual model and then calculate the bounds with the model. However, it is not necessary for the Kriging model to be accurate everywhere in its input space, because it will need more training samples and hence decrease the efficiency. Since the objective is to determine those bounds, we only need the Kriging model to be partially accurate near the regions of interest, as long as it can help find those bounds accurately. This way, we can save more computational efforts. To this end, we develop a task-oriented learning function and a stopping criterion to adaptively train the Kriging model. We start with an analytic model for the one-dimensional (1D) photoelectron sheath near the lunar surface [3,4]. This model is computationally cheap, and hence, the accurate results can be obtained by brute force. With the accurate results, we can test the accuracy of the proposed method. It is noted here that the ultimate application of this method is not the simple 1D problem presented in this work, but more complicated or computationally expensive models such as three-dimensional (3D) fully kinetic particle-in-cell plasma simulations.

The rest of this paper is organized as follows. Section 2 presents the 1D photoelectron sheath and dust levitation problem on lunar surface, as well as the 1D analytic model. Section 3 briefly introduces the Kriging method and general AKM. Section 4 presents the proposed AKM. Section 5 presents the results. Conclusions are given in Sec. 6.

## 2 Problem Statement

### 2.1 One-Dimensional Photoelectron Sheath Model on the Lunar Surface.

We employ the recently derived 1D photoelectron sheath model for the lunar surface [3,4]. As given in detail in Refs. [3] and [4], there are three types of electric potential profiles [3–6] in the photoelectron sheath: type A, type B, and type C, as shown in Fig. 1, where $\varphi $ is the electric potential and $Z$ is the vertical coordinate. In this study, we focus on type C sheath profile as it is expected at the polar regions of the Moon, where the next lunar landing mission will likely occur.

A typical type C sample curve of $E(Z;P)$ is shown in Fig. 2. Note that both $\varphi $ and $E$ converge to zero at large values of $Z$ where it is used as the electric potential reference (zero potential and zero field).

### 2.2 Dust Levitation.

where $r$ is the radius of the lunar dust grain and $\rho =1.8\u2009g/cm3$ is the mass density of dust grains [10].

### 2.3 Objective.

Due to the lack of information, it is almost impossible to obtain the distribution functions of $P$. The bounds of $P$, however, are much easier to obtain. In some work designs on lunar surface, we need to determine the bounds of $\varphi (Z;P)$ and/or $E(Z;P)$, given the bounds of $P$. In this study, all the parameters in $P$ are modeled as interval variables, whose domain is denoted as $\Omega $. For a given realization $p$ of $P$, both $\varphi (Z;p)$ and $E(Z;p),Z\u2208[Zmin,Zmax]$ are obtained by solving the ODE.

However, this method is too expensive or even unaffordable. One reason is that solving the ODE a large number $NMCS$ of times is time-consuming, even when the analytic solution to the ODE is available for the 1D problem. Another reason is that there is no analytic solution to complex two-dimensional or 3D problems where kinetic particle-in-cell simulations are usually employed to solve the electrostatic field through Poisson's equation.

The objective of this study is to develop a method to determine $\varphi \xaf(Z)$, $\varphi \xaf(Z)$, $E\xaf(Z)$, and $E\xaf(Z)$ accurately and then calculate $Z*$ of dust grains. It is noted here that the ultimate application of this method is not the relatively simple 1D problem presented in this work but more complicated or computationally expensive models such as 3D fully kinetic particle-in-cell plasma simulations. For computationally expensive models, evaluating the model consumes the majority of computational resource, so we will use the number $NODE$ of ODEs that we need to solve as a measure of the computational cost.

## 3 Introduction to Kriging Model and Adaptive Kriging Method

### 3.1 Overview of Kriging Model.

Kriging model makes regression to a black-box function (BBF) using a training sample set, or a design of experiment (DoE). The main idea of Kriging is to treat the BBF as a realization of a Gaussian random field indexed by the input variables of the BBF. The theoretical foundation of Kriging model is exactly the Bayesian inference [29]. From the perspective of Bayesian interface, a prior Gaussian random field is trained by the DoE, and hence, a posterior Gaussian random field is generated. Then the mean value function, also indexed by the input variables of the BBF, of the posterior random field is the Kriging prediction to the BBF. Besides, the variance function, also indexed by the input variables of the BBF, of the posterior random field quantifies the local prediction uncertainty or prediction error.

The randomness, or uncertainty, of the posterior random field mainly comes from the fact that only a limit number of samples of the BBF are used to train the prior random field. In other words, only part of the information of the BBF is available, and the missing part of information leads to the epistemic uncertainty in the random field. Generally, the more training samples we use, the less epistemic uncertainty will result and with stronger confidence will we predict the BBF.

### 3.2 Formulation for Kriging Model.

where $\mu $ is an unknown parameter representing the mean value of the random field $K(X)$ and $\eta (X;\xi 2,\theta )$ is a zero-mean stationary Gaussian random field indexed by $X$, the input variables of a BBF $k(X)$. Both the variance parameter $\xi 2$ and correlation parameters $\theta $ of $\eta (X;\xi 2,\theta )$ are unknown. The parameters $\mu $, $\xi 2$, and $\theta $ fully define the prior random field $K(X)$. A DoE, or a training sample set, of $k(X)$ is used to train $K(X)$ and then determine those parameters.

where $D$ is the dimension of $X$, $xd(i)$ is the $dth$ component of $x(i)$, and $\theta d$ is the $dth$ component of $\theta $.

For a BBF $k(X)$, the Kriging model predicts $k(x)$ as $k\u0303(x)$, which is a normal variable whose mean value and variance are $k\u0302(x)$ and $\sigma 2(x)$, respectively. Note that $\sigma 2(x)$ is also termed as mean squared error. Generally, $k\u0302(x)$ is regarded as the deterministic prediction to $k(x)$, since a deterministic prediction is usually needed. $\sigma 2(x)$ measures the prediction uncertainty, or prediction error, and therefore, we can validate a Kriging model simply using $k\u0302(x)$ and $\sigma 2(x)$ without employing traditional validation methods, such as the cross validation [30]. Because of this advantage, many algorithms have been proposed to adaptively train a Kriging model for expensive BBFs [14–27,31–36]. When sufficient training samples have been used for training, $\sigma 2(x)$ converges to 0, and the normal variable $k\u0303(x)$ degenerates to a deterministic value, i.e., the exact value of $k(x)$.

### 3.3 An Example of Kriging Model.

Figure 4 shows a 1D example of Kriging model. In total, five initial training samples are used to train the Kriging. The vertical distance between $k\u0302(x)\xb1\sigma (x)$ graphically quantify the prediction error at $x$. The larger the distance, the larger the prediction error. On interval $[0,\u20092]$, the training samples are denser than that on $[2,\u200910]$. Consequently, the prediction error is smaller on $[0,\u20092]$ than that on $[2,\u200910]$. It is noted that the prediction error is not only dependent on the density of the training samples but also the nonlinearity of the BBF. With the prediction error shown in Fig. 4, it is obvious that in order to improve the prediction accuracy, we need to add training samples somewhere near $x=4$ and $x=8$. Figure 5 shows the updated Kriging model with one more training sample added at $x=8$. The overall prediction accuracy is improved significantly.

### 3.4 Adaptive Kriging Method.

The main idea of AKM is to adaptively add training samples to update the Kriging model iteratively until an expected accuracy is achieved. Figure 6 shows a brief flowchart of AKM. The QoI is what we aim to calculate, such as $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$. Since the QoI is calculated through the Kriging model instead of the BBF itself, there is an inevitable error caused by the Kriging model. The error metric is used to measure the error. The stopping criterion, which is based on the error metric, is used to determine when to stop adding training samples. Once the error of QoI is sufficiently small, it is reasonable to return the QoI and stop the algorithm. If the error is large in an iteration, we must add one or more training samples to update the Kriging model. How to determine new training samples is the task of the learning function. A good learning function should be robust and lead to a high convergence rate.

Given a specific engineering problem, the key of employing an AKM is to make good use of all available information, such as the features of the BBF and QoI, and then design a customized or task-oriented error metric, stopping criterion and learning function.

## 4 The Proposed Method

In this section, we present detailed procedures of calculating $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$. Similar procedures can also apply to calculate $E\xaf(Z)$ and $E\xaf(Z)$.

### 4.1 Overview of the Proposed Method.

In step 5, an error metric is developed to measure the error of $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ estimated by Eqs. (16) and (17). Step 6 is about a stopping criterion. Details about steps 5 and 6 will be given in Subsec. 4.4. The learning function involved in step 7 will be given in Subsec. 4.3. The implementation of the proposed method will be given in Subsec. 4.5.

There are two significant differences between most existing AKMs and the proposed method. First, the former aims at estimating a constant value, such as the structural reliability and robustness, while the latter aims at estimating two functions, i.e., $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$. Second, when given a specific value of input, the output of the BBFs involved in the former methods is a single value. However, in this work, with a given realization $p$ of $P$, the output of solving the ODE is a function $\varphi (Z;p)$. With those differences, we cannot use the existing error metrics, stopping criteria or learning functions. Instead, we take into consideration those differences and design a new error metric, stopping criterion and learning function to fit the problem. This is the main contribution of the proposed algorithm.

### 4.2 Candidate Samples and Initial Training Samples.

For numerical computation, we need to evenly discretize $\Omega $ into a few points. Suppose $Pi$, the *i*$th$ component of $P$, is discretized into $Ni$ points, then $\Omega $ will be discretized into in total $NP=\u220fi=15Ni$ points. For convenience, we denote the set of those points as $pMCS$. Similarly, $Z$ is discretized into $NZ$ points (denoted as $zMCS$) in its range $[Zmin,Zmax]$. Theoretically, any $p\u2208\Omega $ could be selected as a training sample for $\varphi \u0302(Z;P)$. However, we do not want any two training samples to be clustering together, because we use the exact interpolation in Kriging and clustered training samples may impact the training and the convergence rate of the proposed AKM. Therefore, we only select training samples of $P$ from $pMCS$ and call $pMCS$ candidate samples or candidate points.

where $pi(h)$ is the $ith$ component of $p(h)$, $pi(c)$ is the $ith$ component of $p(c)$, $Pi,max$ is the maximal value of $Pi$ which is the $ith$ component of $P$, and $Pi,min$ is the minimal value of $Pi$. Then $p(h)$ is rounded to $p*=argminp\u2208pMCSd(p(h),p(c))$. When all the initial training samples generated by Hammersley have been rounded to the nearest ones in $pMCS$, we get the initial training sample set $pin\u2282pMCS$ of $P$.

Solving the ODE $Nin$ times, each with a sample in $pin$, we get $Nin$ samples of $\varphi (Z;P)$. Note that each sample of $\varphi (Z;P)$ has $NZ$ points, since we discretized $Z$ into $NZ$ points. Then, we have $NZNin$ input training points $zMCS\xd7pin$. Except the $Nin$ points at $Zmax$, we select the other $(NZ\u22121)Nin$ points to form the first part of input training sample set of $\varphi \u0302(Z;P)$. We denote those $(NZ\u22121)Nin$ input training points as $xinp1$, where superscript inp of $x$ represents input, and the superscript *1* means that $xinp1$ is only the first part of the entire input training sample set. The other part $xinp2$ is given below.

The input training sample set $xinp=xinp1\u222axinp2$. Denote the corresponding electric potential $\varphi $ at $xinp$ as $\varphi out$. The input–output training sample pairs $(xinp,\varphi out)$ are used to build the initial $\varphi \u0302(Z;P)$. More training samples will be added to update $\varphi \u0302(Z;P)$ later.

### 4.3 Learning Function.

Generally, the initial Kriging model is not accurate enough to get $\varphi \xaf(Z)$ or $\varphi \xaf(Z)$ accurately through Eqs. (5) and (6). To improve the accuracy of $\varphi \u0302(Z;P)$ and hence of $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$, we need to add training samples of $\varphi (Z;P)$ to refine $\varphi \u0302(Z;P)$. A learning function is used to determine which sample of $P$, and hence of $\varphi (Z;P)$, should be added.

where $p(next)$ is the next to-be-added sample of $P$, $\varphi \u0302(z;p)$ is the predicted value of $\varphi (z;p)$ by the Kriging model $\varphi \u0302(Z;P)$, and $\sigma (z;p)$ is the standard deviation of the prediction. Both $\varphi \u0302(z;p)$ and $\sigma (z;p)$ are calculated by the Kriging toolbox. $\sigma (z(j);p)\varphi \u0302(z(j);p)$ is the deviation coefficient of the prediction at $(z;p)$, and thus, the learning function in Eq. (20) determines the training sample $p(next)$ at which the summation of the absolute deviation coefficients of the predictions along $Z$ coordinate is maximal. The summation $\u2211z\u2208zMCS|\sigma (z;p)\u2009\varphi \u0302(z;p)|$ measures the overall prediction error at $p$. Adding a sample of $\varphi (Z;P)$ at $p$ to update $\varphi \u0302(Z;P)$ will let $\u2211z\u2208zMCS|\sigma (z;p)\u2009\varphi \u0302(z;p)|$ become zero, and therefore adding a sample of $\varphi (Z;P)$ at $p(next)$ to update $\varphi \u0302(Z;P)$ will decrease the overall prediction error of $\varphi (Z;P)$ by the largest extent. This is the basic mechanism of the learning function in Eq. (20).

However, we do not necessarily need $\varphi \u0302(Z;P)$ to be overall accurate. Since the objective is to estimate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ accurately and efficiently, we only need $\varphi \u0302(Z;P)$ to be partially or locally accurate enough so that it can help estimate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ accurately. With this idea, we can further improve the efficiency of updating $\varphi \u0302(Z;P)$ by adding training samples more skillfully.

where we sum up the absolute values of the improvement rate. This learning function means that if we added a training sample $\varphi (Z;p(next))$, which has $NZ$ points, to update $\varphi \u0302(Z;P)$, we could expect to get the best improvement of $\varphi \xaf(Z)$.

Once $p(next)$ has been determined, we solve the ODE to numerically get a function $\varphi (Z;p(next))$, from which we get $(NZ\u22121)$ points (the remaining one at $Zmax$, where $\varphi \u22610$, is excluded) and add them into $(xinp,\varphi out)$ to enrich the training samples.

### 4.4 Error Metric and Stopping Criterion.

where $\gamma $ is a threshold that controls the efficiency and accuracy of the proposed AKM. Generally speaking, a smaller $\gamma $ will lead to higher accuracy but lower efficiency.

### 4.5 Implementation.

where $mean(\xb7)$ represents mean value.

Based on all the procedures given above, we generate the pseudo-codes of the proposed AKM given in Algorithm 1.

Row | Pseudo-codes |
---|---|

1 | Evenly discretize $\Omega $ into $NP$ points $pMCS$. |

2 | Evenly discretize interval $[Zmin,Zmax]$ into $NZ$ points $zMCS$. |

3 | Generate $Nin$ samples $pin$ of $P$ with procedures given in Subsec. 4.2. |

4 | Solve ODE $Nin$ times to get $Nin$ samples $\varphi (Z;p),p\u2208pin$ of $\varphi (Z;P)$; Calculate $\u03f5$ with Eq. (29); $NODE=Nin$. |

5 | Determine $(xinp,\varphi out)$ with procedures given in Subsec. 4.2; $\varphi out=\varphi out+\u03f5$. |

6 | WHILE TRUE DO |

7 | Build Kriging model $\varphi (Z;P)$ using $(xinp,\varphi out)$. |

8 | Calculate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ with Eqs. (16) and (17); $\varphi \xaf(Z)=\varphi \xaf(Z)\u2212\u03f5$; $\varphi \xaf(Z)=\varphi \xaf(Z)\u2212\u03f5$. |

9 | Calculate $\Gamma $ with Eq. (27). |

10 | IF$(\Gamma \u2265\gamma )$DO |

11 | Solve Eq. (20) for $p(next)$; $NODE=NODE+1.$ |

12 | Solve ODE to get a new sample $\varphi (Z;p(next))$; $\varphi (Z;p(next))=\varphi (Z;p(next))+\u03f5$; All points of $\varphi (Z;p(next))$ excluding the one at $Zmax$ are added into $(xinp,\varphi out).$ |

13 | ELSE |

14 | BREAK WHILE |

15 | END IF |

16 | END WHILE |

17 | RETURN$\varphi \xaf(Z)$, $\varphi \xaf(Z)$, and $NODE$. |

Row | Pseudo-codes |
---|---|

1 | Evenly discretize $\Omega $ into $NP$ points $pMCS$. |

2 | Evenly discretize interval $[Zmin,Zmax]$ into $NZ$ points $zMCS$. |

3 | Generate $Nin$ samples $pin$ of $P$ with procedures given in Subsec. 4.2. |

4 | Solve ODE $Nin$ times to get $Nin$ samples $\varphi (Z;p),p\u2208pin$ of $\varphi (Z;P)$; Calculate $\u03f5$ with Eq. (29); $NODE=Nin$. |

5 | Determine $(xinp,\varphi out)$ with procedures given in Subsec. 4.2; $\varphi out=\varphi out+\u03f5$. |

6 | WHILE TRUE DO |

7 | Build Kriging model $\varphi (Z;P)$ using $(xinp,\varphi out)$. |

8 | Calculate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ with Eqs. (16) and (17); $\varphi \xaf(Z)=\varphi \xaf(Z)\u2212\u03f5$; $\varphi \xaf(Z)=\varphi \xaf(Z)\u2212\u03f5$. |

9 | Calculate $\Gamma $ with Eq. (27). |

10 | IF$(\Gamma \u2265\gamma )$DO |

11 | Solve Eq. (20) for $p(next)$; $NODE=NODE+1.$ |

12 | Solve ODE to get a new sample $\varphi (Z;p(next))$; $\varphi (Z;p(next))=\varphi (Z;p(next))+\u03f5$; All points of $\varphi (Z;p(next))$ excluding the one at $Zmax$ are added into $(xinp,\varphi out).$ |

13 | ELSE |

14 | BREAK WHILE |

15 | END IF |

16 | END WHILE |

17 | RETURN$\varphi \xaf(Z)$, $\varphi \xaf(Z)$, and $NODE$. |

### 4.6 Validation Discussion.

Theoretically, it is vital to validate the Kriging model to make sure that it has been trained accurately. An explicit validation, however, is not involved in the proposed AKM. There are two main reasons. First, the adaptive training focuses on the accuracy of QoI instead of the accuracy of the Kriging model. Once there is an indication that the accuracy of QoI in current training iteration is sufficient, i.e., the stopping criterion in Eq. (27) is satisfied, the algorithm stops adding more training samples, no matter the Kriging model is globally accurate or not. As a result, when the algorithm has converged, it is very likely that the Kriging model is accurate only on some subdomains but not accurate on other domains. Therefore, it is not suitable to do explicit cross validation. Second, the error metric $\Gamma $ can measure the accuracy of QoI, and therefore we in fact do validation implicitly. As long as the accuracy of QoI is sufficient, it does not matter if the Kriging model is or not accurate on the entire domain.

## 5 Results

In this section, we illustrate the proposed AKM. MCS is used to solve the same problems with brute force. Results by MCS are treated as standard to verify the proposed AKM. We build the Kriging model and calculate the Kriging predictions using the DACE toolbox [39]. The anisotropic Gaussian kernel is used.

### 5.1 Sheath Profile.

We consider the 1D photoelectron sheath problem discussed in Sec. 2. The sun elevation angle is given as 9 deg. The maximal and minimal values of $P=(vd,Te,Tp,ni,\u221e,np)$ are given in Table 1. We use both MCS and the proposed AKM to estimate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$. The values of all involved parameters are given in Table 2.

Variables | $vd\u2009(m/s)$ | $Te\u2009(eV)$ | $Tp\u2009(eV)$ | $ni,\u221e\u2009(cm\u22123)$ | $np\u2009(cm\u22123)$ |
---|---|---|---|---|---|

Minimum | $421,200$ | $10.8$ | $1.8$ | $7.83$ | $57.6$ |

Maximum | $414,800$ | $13.2$ | $2.2$ | $9.57$ | $70.4$ |

Variables | $vd\u2009(m/s)$ | $Te\u2009(eV)$ | $Tp\u2009(eV)$ | $ni,\u221e\u2009(cm\u22123)$ | $np\u2009(cm\u22123)$ |
---|---|---|---|---|---|

Minimum | $421,200$ | $10.8$ | $1.8$ | $7.83$ | $57.6$ |

Maximum | $414,800$ | $13.2$ | $2.2$ | $9.57$ | $70.4$ |

Parameters | $N1\u223cN5$ | $NP$ | $Nin$ | $NP\u2032$ | $NZ$ | $\gamma $ |
---|---|---|---|---|---|---|

Values | $5$ | $55$ | $5$ | $100$ | $50$ | $0.01$ |

Parameters | $N1\u223cN5$ | $NP$ | $Nin$ | $NP\u2032$ | $NZ$ | $\gamma $ |
---|---|---|---|---|---|---|

Values | $5$ | $55$ | $5$ | $100$ | $50$ | $0.01$ |

The domain $\Omega $ of $P$ is discretized into $NP=55$ points, which are assembled into $pMCS$. The $Nin=5$ samples in hypercube space $[0,1]5$, generated by Hammersley sampling method, are given in Table 3. Then, the five samples are mapped into $\Omega $, as given in Table 4. Rounding the five samples in $\Omega $ to the nearest ones in $pMCS$, we get the initial samples $pin$ of $P$, as given in Table 5. Solving the ODE five times, each with a sample in $pin$, we get five samples of $\varphi (Z;P)$ as shown in Fig. 8.

Sample number | Dimension 1 | Dimension 2 | Dimension 3 | Dimension 4 | Dimension 5 |
---|---|---|---|---|---|

1 | 0 | 0.5000 | 0.3333 | 0.2000 | 0.1429 |

2 | 0.2 | 0.2500 | 0.6667 | 0.4000 | 0.2857 |

3 | 0.4 | 0.7500 | 0.1111 | 0.6000 | 0.4286 |

4 | 0.6 | 0.1250 | 0.4444 | 0.8000 | 0.5714 |

5 | 0.8 | 0.6250 | 0.7778 | 0.0400 | 0.7143 |

Sample number | Dimension 1 | Dimension 2 | Dimension 3 | Dimension 4 | Dimension 5 |
---|---|---|---|---|---|

1 | 0 | 0.5000 | 0.3333 | 0.2000 | 0.1429 |

2 | 0.2 | 0.2500 | 0.6667 | 0.4000 | 0.2857 |

3 | 0.4 | 0.7500 | 0.1111 | 0.6000 | 0.4286 |

4 | 0.6 | 0.1250 | 0.4444 | 0.8000 | 0.5714 |

5 | 0.8 | 0.6250 | 0.7778 | 0.0400 | 0.7143 |

Sample number | $vd\u2009(m/s)$ | $Te\u2009(eV)$ | $Tp\u2009(eV)$ | $ni,\u221e\u2009(\u2009cm\u22123)$ | $np\u2009(\u2009cm\u22123)$ |
---|---|---|---|---|---|

1 | 421,200 | 12.0000 | 1.9333 | 8.1780 | 59.4286 |

2 | 439,920 | 11.4000 | 2.0667 | 8.5260 | 61.2571 |

3 | 458,640 | 12.6000 | 1.8444 | 8.8740 | 63.0857 |

4 | 477,360 | 11.1000 | 1.9778 | 9.2220 | 64.9143 |

5 | 496,080 | 12.3000 | 2.1111 | 7.8996 | 66.7429 |

Sample number | $vd\u2009(m/s)$ | $Te\u2009(eV)$ | $Tp\u2009(eV)$ | $ni,\u221e\u2009(\u2009cm\u22123)$ | $np\u2009(\u2009cm\u22123)$ |
---|---|---|---|---|---|

1 | 421,200 | 12.0000 | 1.9333 | 8.1780 | 59.4286 |

2 | 439,920 | 11.4000 | 2.0667 | 8.5260 | 61.2571 |

3 | 458,640 | 12.6000 | 1.8444 | 8.8740 | 63.0857 |

4 | 477,360 | 11.1000 | 1.9778 | 9.2220 | 64.9143 |

5 | 496,080 | 12.3000 | 2.1111 | 7.8996 | 66.7429 |

Sample number | $vd\u2009(m/s)$ | $Te\u2009(eV)$ | $Tp\u2009(eV)$ | $ni,\u221e\u2009(\u2009cm\u22123)$ | $np\u2009(\u2009cm\u22123)$ |
---|---|---|---|---|---|

1 | 421,200 | 12.0000 | 1.9000 | 8.2650 | 60.8000 |

2 | 444,600 | 11.4000 | 2.1000 | 8.7000 | 60.8000 |

3 | 468,000 | 12.6000 | 1.8000 | 8.7000 | 64.0000 |

4 | 468,000 | 11.4000 | 2.0000 | 9.1350 | 64.0000 |

5 | 491,400 | 12.0000 | 2.1000 | 7.8300 | 67.2000 |

Sample number | $vd\u2009(m/s)$ | $Te\u2009(eV)$ | $Tp\u2009(eV)$ | $ni,\u221e\u2009(\u2009cm\u22123)$ | $np\u2009(\u2009cm\u22123)$ |
---|---|---|---|---|---|

1 | 421,200 | 12.0000 | 1.9000 | 8.2650 | 60.8000 |

2 | 444,600 | 11.4000 | 2.1000 | 8.7000 | 60.8000 |

3 | 468,000 | 12.6000 | 1.8000 | 8.7000 | 64.0000 |

4 | 468,000 | 11.4000 | 2.0000 | 9.1350 | 64.0000 |

5 | 491,400 | 12.0000 | 2.1000 | 7.8300 | 67.2000 |

Each sample of $\varphi (Z;P)$ contains $NZ=50$ numerical points. Excluding the five points at $Zmax$, we have $NZNin\u22125=245$ training points in $(\u2009xinp1,\varphi out1)$. With Hammersley sampling method, we generate $NP\u2032=100$ samples of $P$ and hence $100$ training points in $(xinp2,\varphi out2)$. Note that all points in $(xinp2,\varphi out2)$ have $Z=Zmax$ and $\varphi =0$. Combining $(\u2009xinp1,\varphi out1)$ and $(xinp2,\varphi out2)$, we have 345 training points in $(xinp,\varphi out)$. To do the translation mentioned in Subsec. 4.5, we update $\varphi out$ simply by $\varphi out=\varphi out+\u03f5$, where $\u03f5\u2009=\u2009\u22126.97\u2009V$ is obtained with Eq. (29). With the updated $(xinp,\varphi out)$, we build an initial Kriging model and then estimate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ through Eqs. (16) and (17). Finally, we translate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ back by $\varphi \xaf(Z)=\varphi \xaf(Z)\u2212\u03f5$ and $\varphi \xaf(Z)=\varphi \xaf(Z)\u2212\u03f5$. Figure 9 shows the $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ estimated by both MCS and the proposed AKM (with the initial Kriging model). It shows that the initial Kriging model is not able to predict $\varphi \xaf(Z)$ or $\varphi \xaf(Z)$ with sufficient accuracy.

To improve the accuracy, the proposed method indicates adding a sample at $p(next)=(514800,13.2,\u20092.2,9.57,\u200957.6)$. With the $p(next)$, we solve the ODE and get a new sample of $\varphi (Z;P)$. This sample contains $NZ=50$ numerical points. We translate all the numerical points and add them, excluding the one at $Zmax$, to update $(xinp,\varphi out)$. The reason why we abandon the point at $Zmax$ is that there are already enough points at $Zmax$ in $(xinp2,\varphi out2)$. With the updated $(xinp,\varphi out)$, we build a new $\varphi \u0302(Z;P)$. With the new $\varphi \u0302(Z;P)$ another $p(next)$ is indicated. With similar procedures, more and more samples of $\varphi (Z;P)$ are added to refine $\varphi \u0302(Z;P)$ until the stopping criterion given in Eq. (28) is satisfied.

The final estimation of $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ is shown in Fig. 10. It shows that the proposed AKM can estimate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$ very accurately. 16 more samples of $\varphi (Z;P)$ have been added to refine $\varphi \u0302(Z;P)$, and therefore, in total $NODE=Nin+16=21$ ODE solutions are needed. Compared to $NP=3125$ ODE solutions needed in MCS, the proposed method is very efficient.

### 5.2 Dust Levitation.

In this example, we still consider the same 1D photoelectron sheath problem in Subsec. 5.1, but the objective is to estimate $E\xaf(Z)$ and $E\xaf(Z)$ and then calculate the dust levitation height. The values of all involved parameters are given in Table 6.

Parameters | $N1\u223cN5$ | $NP$ | $Nin$ | $NP\u2032$ | $NZ$ | $\gamma $ |
---|---|---|---|---|---|---|

Values | $5$ | $55$ | $5$ | $100$ | $50$ | $0.01$ |

Parameters | $N1\u223cN5$ | $NP$ | $Nin$ | $NP\u2032$ | $NZ$ | $\gamma $ |
---|---|---|---|---|---|---|

Values | $5$ | $55$ | $5$ | $100$ | $50$ | $0.01$ |

The procedures used to estimate $E\xaf(Z)$ and $E\xaf(Z)$ are almost the same as that used to estimate $\varphi \xaf(Z)$ and $\varphi \xaf(Z)$. The only difference is that the samples of $E(Z;P)$ instead of $\varphi (Z;P)$ are used. The final estimation of $E\xaf(Z)$ and $E\xaf(Z)$ is shown in Fig. 11. It shows that the proposed AKM method is very accurate. As for the efficiency, the proposed method needs only $NODE=Nin+18=23$ ODE solutions. Compared to $NP=3125$ ODE solutions needed in MCS, the proposed method is very efficient.

When the upper and lower bounds of the electric field have been determined, we can use them to determine the levitation heights of the dust grains. Assuming there are two types of dust grains, A and B, in the electric field. The relevant parameters of the grains are given in Table 7, where $e=1.062\xd710\u221219\u2009C$ is the electric charge of an electron. The dust levitation heights are shown in Fig. 12 and given in Table 8. Due to the uncertainty of $P$, the levitation heights of both A and B are also uncertain. The levitation height of A may be any value in the interval $[0\u2009m,\u20099.33\u2009m]$, which is estimated by the proposed method. The interval determined by MCS is $[0\u2009m,\u20099.26\u2009m]$. It shows that the proposed method can estimate the levitation height of grain A with sufficient accuracy. Similar conclusion applies to the levitation height of grain B.

Grains | $r\u2009(\mu m)$ | $m\u2009(g)$ | $q/e$ | $w\u2009(V/m)$ |
---|---|---|---|---|

A | $0.5$ | $1.5268\xd710\u221212$ | $50,000$ | $\u22120.4658$ |

B | $0.3$ | $3.2979\xd710\u221213$ | $45,000$ | $\u22120.1118$ |

Grains | $r\u2009(\mu m)$ | $m\u2009(g)$ | $q/e$ | $w\u2009(V/m)$ |
---|---|---|---|---|

A | $0.5$ | $1.5268\xd710\u221212$ | $50,000$ | $\u22120.4658$ |

B | $0.3$ | $3.2979\xd710\u221213$ | $45,000$ | $\u22120.1118$ |

Grains | AKM | MCS | Relative error (%) | |
---|---|---|---|---|

A | $Zmin*\u2009(m)$ | $0.00$ | 0.00 | $0.0$ |

$Zmax*\u2009(m)$ | $9.33$ | 9.26 | $0.8$ | |

B | $Zmin*\u2009(m)$ | $10.88$ | 11.00 | $\u22121.1$ |

$Zmax*\u2009(m)$ | $25.55$ | 25.55 | $0.0$ |

Grains | AKM | MCS | Relative error (%) | |
---|---|---|---|---|

A | $Zmin*\u2009(m)$ | $0.00$ | 0.00 | $0.0$ |

$Zmax*\u2009(m)$ | $9.33$ | 9.26 | $0.8$ | |

B | $Zmin*\u2009(m)$ | $10.88$ | 11.00 | $\u22121.1$ |

$Zmax*\u2009(m)$ | $25.55$ | 25.55 | $0.0$ |

Given any dust grain with known $w$ value, we can easily determine its levitation height interval using the method shown in Fig. 12. This will help to evaluate the risk or damage caused by the levitated dust grains for lunar exploration missions.

## 6 Conclusions

We presented an adaptive Kriging-based method to perform UQ analysis of the 1D photoelectron sheath and dust levitation on the lunar surface. A recently derived 1D photoelectron sheath model was used as the high-fidelity physics-based model and the black-box function. Adaptive Kriging method, with a task-oriented learning function and stopping criterion, was utilized to improve the efficiency in calculating the upper and lower bounds of electric potential as well as dust levitation height, given the intervals of model parameters. Experiment analysis shows that the proposed AKM method is both accurate and efficient. Current and ongoing efforts are focused on building adaptive Kriging model for two-dimensional and 3D kinetic particle simulations of lunar plasma/dust environment and perform UQ analysis.

## Acknowledgment

We would like to acknowledge the support from NASA-Missouri Space Grant Consortium through NASA-EPSCoR Missouri and NSF-CMMI No. 1923799.

## Funding Data

National Aeronautics and Space Administration (NASA-Missouri Space Grant Consortium) (Funder ID: 10.13039/100000104).

National Science Foundation (CMMI No. 1923799; Funder ID: 10.13039/100000001).

## References

**59**(1), pp.

**54**(3), pp.

**10**(3), pp.

**61**(3), pp.