{'fraction_non_alphanumeric': 0.03255829358770535, 'fraction_numerical': 0.015103338632750398, 'mean_word_length': 4.485646802325581, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simplification of Carbon Bond Mechanism IV (CBM-IV) under Different Initial Conditions by Using Concentration Sensitivity Analysis

Carbon Bond Mechanism IV (CBM-IV) is a widely used reaction mechanism in which VOCs are grouped according to the molecular structure. In the present study, we applied a sensitivity analysis on the CBM-IV mechanism to clarify the importance of each reaction under two different initial conditions (urban and low-NOx scenarios). The reactions that exert minor influence on the reaction system are then screened out from the mechanism, so that a reduced version of the CBM-IV mechanism under specific initial conditions can be obtained. We found that in a typical urban condition, 11 reactions can be removed from the original CBM-IV mechanism, and the deviation is less than 5% between the results using the original CBM-IV mechanism and the reduced mechanism. Moreover, in a low-NOx initial condition, two more reactions, both of which are nitrogen-associated reactions, can be screened out from the reaction mechanism, while the accuracy of the simulation is still maintained. It is estimated that the reduction of the CBM-IV mechanism can save 11–14% of the computing time in the calculation of the chemistry in a box model simulation.


Introduction
Air pollution is becoming more severe at many parts of the world, mainly due to increasing human activities such as the burning of fossil fuels [1]. Many types of pollutants such as carbon monoxide (CO) and sulfur dioxide (SO 2 ) as well as small particles of soot are formed from the incomplete burning of hydrocarbons in power plants or factories and then released into the ambient air, leading to poor air quality in cities and a damage to the health of the people. For example, it was reported by Chan and Yao [2] that the air quality in mega cities of China exceeds the Chinese Grade-II standard on 10-30% of days based on a few years of data, and the occurrence of ozone pollution events also becomes more frequent due to the increased emission of VOCs. It was also reported by the World Health Organization (WHO) that the annual mean levels of the particulate matter (PM) have increased by more than 8% from 2010 to 2015 based on the data collected in approximately 800 cities around the world [3]. Moreover, it was estimated by WHO [3] that in the year 2012, about 3 million deaths were attributed to air pollution from particulate matter. It is now well known that these pollution events are usually caused by an occurrence of a series of chemical reactions in the atmosphere. For example, the photochemical smog occurred in the troposphere arises from an excessive ozone (O 3 ) generated by the oxidation of VOCs in the presence of OH. Thus, in order to understand the principles of the pollution events, formulate control strategies, and solve environmental problems, an accurate description of the complex atmospheric chemical processes is crucial. Atmospheric reaction mechanism, which is an essential part of the air quality model, is a means of representing Therefore, in the present study, we applied a concentration sensitivity analysis on the CBM-IV mechanism to identify important reactions in the mechanism under different initial conditions (typical urban and low-NO x scenarios). Then, based on the results of the concentration sensitivity analysis, we removed the reactions that were indicated unimportant for the change of all the chemical species. By doing that, a reduced version of the CBM-IV mechanism can be obtained, which is able to save the computing time under the specific initial conditions. Moreover, the accuracy of the original CBM-IV mechanism can be kept.
The structure of the manuscript is as following. In the next section (Section 2), the method for the simplification of the CBM-IV mechanism and the configurations of the computations are given. In Section 3, most of the computational results are shown, and a comparison between the results of the simulations by using the original CBM-IV mechanism and the reduced one is made. In the last section (Section 4), major conclusions of this paper are drawn, and an extension of the present study is discussed.

Mathematical Models and Methods
The strategy of the present study is as following. At first, a box model KINAL [29] was used to capture the temporal evolution of the chemical species included in CBM-IV under different initial conditions. The obtained mixing-ratio curves over time were also compared with that obtained by using another open-sourced software KPP [18,19] for the purpose of validation. After verifying the correctness of the KINAL simulations, we continued to perform a concentration sensitivity analysis on the CBM-IV mechanism, followed by a mechanism simplification. By doing that, we were able to obtain a reduced version of the CBM-IV mechanism with less reactions compared to the original one. These procedures were also made by using the box model KINAL. At last, after the simplification of the mechanism, the reduced version of the CBM-IV mechanism was again implemented in KINAL so that a comparison between the results achieved before and after the simplification can be made.
In the following, the details of the model used in the present study as well as the configurations of the computations are given.

Reaction System
In the present study, the change of each chemical species in the reaction system can be described as: In Equation (1), c denotes a column vector of the species concentrations, and c| t=0 = c 0 stands for the initial conditions. The variables, k, t and E, in Equation (1) represent a vector of reaction rate constants, time, and the local surface emissions, respectively. By solving Equation (1), we were able to capture the change of each species over time. As the reaction rates in the CBM-IV mechanism differ a lot, the Jacobi matrix formed during the process of solving Equation (1) is stiff which adds the difficulty in finding a trade-off between the efficiency and the accuracy. In this study, Equation (1) was solved in a box model KINAL [29]. KINAL is an open-sourced program package written in FORTRAN language, aiming for the simulation of chemical kinetic mechanisms. A subroutine of KINAL, DIFF, can solve stiff equations so that the temporal change of the reaction system can be derived.
After solving Equation (1), we validated the obtained results by using another open-sourced software KPP [18,19]. KPP is a chemical box model for the investigation of dynamic chemical systems, which has been successfully applied in the study of the tropospheric and stratospheric chemistry [18,30,31]. The reason we chose KPP for validation is that many commonly used reaction mechanisms such as CBM-IV and SAPRC-99 [10,11] are originally included in KPP. In our study, we found the deviation between the results of these two different models (KINAL and KPP) with same configurations lower than 1% (see the supplementary material of this manuscript and also the reference [32]), which ensures the correctness of the KINAL computations and enables us for a further concentration sensitivity analysis.

Concentration Sensitivity Analysis
After validating the KINAL results, we continued to conduct a concentration sensitivity analysis on CBM-IV to reduce the size of the mechanism under different initial conditions. In the concentration sensitivity analysis, the importance of the j-th reaction for the i-th component can be identified by the relative concentration sensitivity coefficient S ij . It is expressed as: In Equation (2), c i , k j represent the concentration of the i-th constituent, and the rate constant of the j-th reaction, respectively. S ij = ∂c i /∂k j denotes the absolute concentration sensitivity, of which the unit depends on the order of the j-th reaction. To compare the value of sensitivity coefficients belonging to different reactions, the normalized sensitivity coefficient, S ij , is introduced by multiplying S ij with c i /k j , so that S ij is a non-dimensional variable. The relative concentration sensitivity S ij can indicate the relative change of the i-th species concentration when a small perturbation occurs on the rate constant of the j-th reaction.
In this study, S ij can be calculated by taking the partial derivative of the i-th constituent in Equation (1) over k j . When the local surface emission E in Equation (1) is assumed independent on the rate constant, Equation (1) becomes: of which the upper limit n s signifies the total number of chemical species included in the mechanism. By replacing ∂c i /∂k j and ∂c l /∂k j in Equation (3) with the absolute concentration sensitivity S ij and S lj , another form of Equation (3) can be obtained: The second term on the right-hand side of Equation (4) denotes the direct change of the i-th species concentration caused by the change of the j-th reaction rate. In contrast to that, the first term on the right-hand side of Equation (4) indicates the indirect effect brought about by the change of the j-th reaction rate. After solving Equation (4), the relative sensitivity can be derived by multiplying the absolute concentration sensitivity S ij with k j /c i . The subroutine SENS in KINAL can calculate a matrix of the absolute concentration sensitivities of the reaction system, which consequently leads to the estimation of the relative sensitivity.
The concentration sensitivity analysis is effective in clarifying the relative importance of a single reaction for a specific component. Accordingly, removing unimportant reactions from the original mechanism can be made. In this research, if the j-th reaction fulfills the criterion shown in Equation (5) according to [33], it is then marked as unimportant and thus can be removed.
n s in Equation (5) is the total number of the chemical components involved, and n t is the total number of time steps in the process of calculation. By using the criterion shown in Equation (5), if the absolute value of the sensitivity belonging to the j-th reaction for all the components is less than 10% at every time point, the j-th reaction can be regarded as unimportant and thus can be removed from the original mechanism.

Configurations of the Model
In this study, we investigated the simplification of the CBM-IV mechanism under two different initial conditions (urban and low-NO x scenarios), and the initial conditions of the model are shown in Table 1. In the simulation under a typical urban condition, the initial mixing ratio of NO x (NO+NO 2 ) is in the order of tens of ppb (ppb = parts per billion). In contrast to that, in the low-NO x case, we reduced the initial NO x amount to 100 ppt (ppt = parts per trillion). Besides, we increased the initial value of HONO from 1 ppb to 30 ppb according to [34]. In the present study, we computed the temporal change of the reaction system within 5 days, and the start time of the simulation is set to 12:00 of the first day.
The reaction mechanism used in the present study is the original version of CBM-IV, which was proposed by Gery et al. [6]. It consists of 33 chemical species and 81 reactions, which are listed in the Appendix A of this manuscript. Generally, four different types of species are included in the CBM-IV mechanism: (1) Inorganic species; (2) Organic species that are explicitly represented due to their unique chemical natures in the environment, such as the formaldehyde (FORM); (3) Organic species denoted by carbon bonds, such as olefins (OLE); (4) Organic species represented by the molecular structure, such as aromatic hydrocarbons (e.g., TOL). The values of the rate constant k in the CBM-IV mechanism used in the present simulation are adopted from Gery et al. [6], and can be found in the Appendix A of this manuscript. They are calculated by using the Arrhenius formula [35]: In Equation (6), R (unit: J K −1 mol −1 ) is the universal gas constant, and E a (unit: J mol −1 ) is the activation energy. T (unit: K) is the temperature, and A represents the pre-exponential factor that has a same unit as k, depending on the order of the reaction. In this study a constant temperature T = 288.15 K is assumed, which is the same as the settings in the KPP examples [34]. Table 1. Initial conditions used in the simulations of the urban scenario [34] and the low-NO x scenario (unit: ppb, ppb = parts per billion). The species not listed in this table are assumed to have a near-zero initial value (10 −5 ppt). With respect to the photolytic reactions in the mechanism, according to KPP [18,19], a radian-like parameter T tmp is introduced to consider the influence on the photolysis rates exerted by the change of the solar zenith angle,

Species
In Equation (7), the sunrise time, T sunrise , is set to 4:30 a.m. in our simulations, and the sunset time T sunset is defined as 7:30 p.m. T local in Equation (7) represents the local hours of the day. When T local resides between 4.5 (the sunrise time) and 19.5 (the sunset time), the photolytic reactions are switched on in the model. Otherwise, they are switched off.

Results and Discussion
To investigate the simplification of the CBM-IV mechanism under different initial conditions, two simulation scenarios were set up, urban and low-NO x cases. The initial condition of the urban scenario was adopted from the benchmark example [34] originally included in KPP (see Table 1). In contrast to that, in the low-NO x scenario, the initial mixing ratios of NO and NO 2 were reduced to 100 ppt, much lower than the typical value range (i.e., tens of ppb) used in the urban simulation.
The simulation results under these two scenarios are presented separately below.

Urban Scenario
In this simulation, the initial air composition contains high levels of pollutants, such as NO x , PAN and O 3 (see Table 1). Figure 1 shows the time behavior of major components, O 3 , NO and NO 2 , under this typical urban condition. These mixing-ratio results over time have been discussed in detail in [32]. Thus, in this paper we only describe them briefly. It is seen in Figure 1a that during the simulated 5 days, the mixing ratio of ozone increases rapidly within the first few hours, from its initial value 100 ppb to a peak value of approximate 180 ppb. After reaching the maximum, the O 3 level starts to decline due to the photolytic decomposition. When the nighttime comes at 7:30 p.m. of the first day (see Figure 1a), the decrease rate of the O 3 mixing ratio becomes lower. However, when the sun rises at 4:30 a.m. of the next day, the amount of ozone continues to drop, until the ozone value reaches a relatively stable level when the end of the simulation is approached.  With respect to the temporal evolution of the nitrogen oxides (NO x = NO + NO 2 ), it is seen in Figure 1a that the variation of NO and NO 2 occurs mostly within the first few hours. Thus, we show the change of the NO x mixing ratios during the first 36 h in Figure 1b. It is seen that the mixing ratio of NO decreases monotonously from its initial value to a near-zero value within 4 h (from 12:00 to 16:00 on the first day). In contrast to the NO profile, the mixing ratio of NO 2 increases at the beginning of the simulation, reaching a peak value of approximately 45 ppb. Then it decreases to less than 1 ppb before hour 16. It should be noted that during the increase of NO 2 , a negative correlation between the mixing ratios of O 3 and NO 2 is found. It is because that at the early stage of the simulation, the amount of NO is abundant. The chemical system composed of ozone and nitrogen oxides is thus mainly controlled by the titration reaction: leading to an enhancement of NO 2 and a decrease of NO and O 3 . However, because the O 3 depleted in Reaction (R1) can be compensated by that formed in the OH radical chain reaction of VOCs: the loss of O 3 is thus found much less than the increase of NO 2 (see Figure 1b). Please note that O( 3 P) in Reaction (R2) denotes the oxygen atom in the ground state. When the reaction proceeds, after hour 16, the amount of NO x given in the initial condition is almost completely consumed due to the conversion of NO 2 to HNO 3 and PAN: As a result, in this time period, the ozone loss and formation is strongly influenced by the photolytic reaction in which O( 1 D) denotes electronically excited state oxygen atom, and the NO x regenerated by the reaction of OH radical with HNO 3 during the daytime: Later, we plotted the values of the ozone sensitivity coefficient corresponding to each reaction in the CBM-IV mechanism, at the time points representing the beginning (the 1st hour) and the end of the simulation (the 96th hour, i.e., 12:00 of the 5th day) (see Figure 2). It is seen in Figure 2a that at the beginning of the simulation, the values of the ozone sensitivity for all the reactions in the mechanism are mostly less than 0.2. Among these reactions, the most dominant reactions at this time are Reactions (SR1) NO 2 + hν → O( 3 P) + NO, (SR3) O 3 + NO → NO 2 and (SR26) OH + NO 2 M − → HNO 3 (see the Appendix A for the index of the reactions). Among these three reactions, it is not surprising that Reactions (SR1) and (SR3) have large influence on the change of the ozone concentration, as these two reactions are parts of the reaction cycle between the nitrogen oxides and ozone. Aside from this, the ozone mixing ratio at this time has the strongest negative dependence on Reaction (SR26), as this reaction is able to convert OH and NO 2 to HNO 3 , leading to a major depletion of ozone during this time period as mentioned above. These findings again suggest that at this time stage, the ozone mixing ratio is mostly determined by the titration reaction and the high initial values of the nitrogen oxides. When the end of the simulation is approached, the ozone sensitivity to each reaction in the mechanism changes (see Figure 2b). In general, most of the O 3 sensitivities increase as the simulation proceeds, denoting an enhanced dependence of the ozone mixing ratio on the local chemistry. It should be noted that Reactions (SR9) and (SR11) have a negative influence on the production of O 3 at the 96th hour (see Figure 2b), while at the first hour these sensitivity coefficients are positive. It is because that at the beginning of the simulation when the amount of NO x is abundant, the photolysis of ozone leads to a formation of OH and a following OH radical chain reaction of VOCs as mentioned above. As a result, the photolytic decomposition of ozone promotes the formation of ozone. In contrast to that, when the end of the 5-day simulation comes, as the nitrogen oxides are almost completely consumed, the amount of ozone formed from the OH radical chain reaction is less than the loss of ozone caused by the photolysis of ozone, leading to a change in the sign of the sensitivities corresponding to these two reactions. Moreover, because the initial nitrogen oxides are mostly converted to HNO 3 and PAN at this end time, ozone is mainly consumed by its photolytic reaction, i.e., Reaction (SR9), instead of the HNO 3 formation reaction (SR26). Thus, the dominance of Reaction (SR9) is highlighted during the end of the simulation (see Figure 2b). In addition, we also found the sensitivity value of Reaction (SR38) declines remarkably over time. This is because that the most deterministic factor for the ozone change at this time is the availability of NO x rather than HO 2 . We then applied the selection criterion described in Equation (5) on the obtained relative concentration sensitivities. Eleven reactions, (SR5), (SR6), (SR20), (SR21), (SR25), (SR40), (SR42), (SR55), (SR56), (SR60) and (SR75) were identified as unimportant so that they can be eliminated from the original reaction mechanism under this urban condition. As a result, a simplified version of the CBM-IV mechanism, consisting of 33 species and 70 reactions, is obtained. We found the change of all the constituents in the simulation using the simplified mechanism is almost identical to that using the original CBM-IV mechanism (not shown here). By comparing the obtained values (see Table 2), we found the maximum deviation between these two results smaller than 5%, which proves the correctness of our mechanism simplification under this urban initial condition. We then estimated the computing time saved by applying the simplified mechanism in KINAL, and it was found that approximately 11% of the computing time can be saved, which enables a faster calculation of the chemistry in applications. However, it should be noted that the exact computing time that can be saved depends on the situation. In our box model computation, the time saved by the mechanism simplification is mainly in the order of minutes to hours, while in 3-D simulations it might be in the order of days or weeks, depending on the mesh resolution and the length of the time step.

Low-NO x Scenario
In this case, we reduced the initial mixing ratio of the nitrogen oxides (NO from 50 to 0.1 ppb. NO 2 from 20 to 0.1 ppb), and increased the initial value of HONO to 30 ppb according to [34]. The temporal evolution of ozone, NO and NO 2 under this low-NO x condition is shown in Figure 3. As seen in Figure 3, the profile change of ozone under these conditions is similar to that in the urban scenario shown in Figure 1. The ozone mixing ratio increases to a peak value within a few hours (see Figure 3a), and then drops to a relatively stable level at the end of the computation. However, the maximum value (∼140 ppb) and the value at the end of the simulation (∼70 ppb) are much lower than those in the urban scenario simulation, 180 ppb and 105 ppb, respectively. These lower ozone values are expected as the formation of ozone is inhibited due to the lack of nitrogen oxides in this low-NO x condition. Moreover, it may also be observed in Figure 3b that under these experimental conditions, the mixing ratio of ozone increases directly from its original value to the peak value, while in the urban scenario the ozone mixing ratio decreases at the beginning (see Figure 1b). The reason for the observed difference may also be attributed to the lower amount of NO and NO 2 under given experimental conditions so that the formation of O 3 by the oxidation of VOCs exceeds the loss of O 3 caused by the titration reaction at the start of the simulation.  Regarding to the nitrogen oxides, it was found that the trends of the NO and NO 2 change are similar in this low-NO x scenario. Both of these two species mixing ratios increase from the start of the computation, reaching a peak value (2.9 ppb for NO 2 and 0.3 ppb for NO) at the first hour (see Figure 3b), and then decline to an amount less than 0.1 ppb at about 2 h after the start of the simulation. The maximum value of the NO 2 mixing ratio in this condition is found remarkably lower than that in the urban scenario, due to the smaller initial value of the nitrogen oxides. Furthermore, different from the urban simulation in which the NO mixing ratio decreases monotonously, in this low-NO x scenario, the value of NO increases at the beginning, due to the photolytic decomposition of the initial HONO at this moment.
The O 3 sensitivity values at the first hour and the 96th hour under low-NO x experimental conditions are shown in Figure 4. Similar to the results of the urban case scenario, at the beginning of the simulation, the O 3 sensitivity is generally lower than 0.2 (see Figure 4a). By comparing Figure 4a with Figure 2a, an enhanced importance of Reaction (SR28) HO 2 + NO → OH + NO 2 is indicated in the low-NO x computation. It is because that when a low value is initially given to the nitrogen oxides, the availability of NO 2 turns out to be the rate-determining factor for the ozone formation. As a result, the ozone becomes very sensitive to the rate change of Reaction (SR28) that converts NO to NO 2 . In contrast to that, at a later period of this simulation (see Figure 4b), the most important reaction for the formation of O 3 is Reaction (SR14) NO 3 + hν → 0.89NO 2 + 0.89O( 3 P) + 0.11NO. It is because that at this time, the NO 3 regenerated from Reaction (SR27) OH + HNO 3 M − → NO 3 during the daytime acts as the major source of the nitrogen oxides as discussed above. In comparison with that, the loss of ozone at this time is strongly influenced by Reactions (SR9), (SR11) and (SR26) (see Figure 4b), which is similar to the situation in the urban simulation shown in Figure 2b.
(a)  In contrast to the urban scenario, in this low-NO x simulation, apart from the 11 reactions identified as negligible in the urban scenario, two more reactions in the mechanism are indicated as unimportant: It was found that during the computation process, the influence of these two reactions on the change of all the constituents in the mechanism is minor, so that they can be eliminated from the mechanism. It was also noticed that both reactions are nitrogen-associated reactions, reflecting a lower importance of the nitrogen related species under this low-NO x initial condition. After removing these 13 reactions indicated in the sensitivity analysis, a reduced reaction mechanism that is composed of 33 species and 68 reactions is obtained under the low-NO x condition. It was observed that the change of all the air components over time by using different reaction mechanisms is similar, denoting a minor difference between the results before and after the simplification procedure. In consistency with the results shown in the previous section, we also see in Table 3 that the deviation of each species mixing ratio is small, with a maximum lower than 3%, which also suggests a successful mechanism reduction. Besides, due to the simplification, it was estimated that 14% of the computing time can be saved in the calculation of the chemistry in this condition.

Conclusions and Future Work
In this study, we used a concentration sensitivity analysis to study the relative importance of individual reaction in the CBM-IV mechanism under different initial conditions. Reactions that exert negligible impact on the reaction system were eliminated from the mechanism so that a mechanism reduction under specific initial conditions can be made. In the present study, we found that when a typical urban condition is initially given, at the beginning of the simulation, the reaction system especially the ozone mixing ratio is deeply influenced by the titration reaction and the high initial values of nitrogen oxides. However, when the end of the 5-day computation is approached, because the nitrogen oxides are mostly converted to HNO 3 and PAN, the ozone value in the CBM-IV mechanism depends heavily on its photolytic reaction and the amount of NO x regenerated by the reaction between HNO 3 and OH. In this condition, we found that 11 reactions are indicated as unimportant and thus can be removed from the original CBM-IV mechanism. The maximum deviation of all the simulated species concentration between the results before and after the simplification was found less than 5%. It was also estimated that the reduction of the mechanism enables a 11% saving of the computing time.
On the contrary, when the initial value of the nitrogen oxides is reduced to 0.1 ppb, the ozone change at an early period of the simulation was found mostly determined by the availability of NO 2 as well as the reaction converting NO to NO 2 . However, when the end of the low-NO x simulation comes, the production of O 3 depends mostly on the photolysis of NO 3 , which is similar to the situation in the urban case. It was also found that 13 reactions including two more nitrogen-associated reactions were identified as negligible in this low-NO x scenario compared to the urban scenario, and can be screened out from the original CBM-IV mechanism. By implementing the mechanism after the simplification into the box model KINAL, we were able to save 14% of the computing time, while the deviation between the results before and after the mechanism reduction is lower than 3%.
In the future, we plan to continue our research by investigating more advanced reaction mechanisms such as CB05 [8], CB6 [9] and SAPRC-07 [12], which are used more frequently at present. The full plan is that we first investigate the internal properties of these reaction mechanisms by using the concentration sensitivity analysis. After that, we will make a simplification of these mechanisms based on the results of the sensitivity analysis and analyze the differences between these mechanisms. Observational data obtained from field campaigns or environmental monitoring stations are also needed to evaluate the mechanisms and confirm the findings. Aside from this, it might also be interesting to figure out the common features of the removed reactions under the given initial conditions. Moreover, the impact brought about by the inclusion of the surface emission on the simplification of the mechanism should be studied. It is also useful to extend the method presented in this manuscript to 3-D simulations. For this purpose, we need to compare the timescale of the atmospheric chemistry with that of air transport such as the horizontal advection and the vertical turbulent mixing. At present, the authors are developing a mechanism simplification method in which the influence of air diffusion can be considered.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1420-3049/24/13/2463/ s1, Table S1: Peak values of major components (NO 2 , H 2 O 2 , O 3 , CO, NO and PAN) obtained in KPP and KINAL, and the deviation of the peak values between these two models. Figure S1: Temporal change of (a) NO 2 , H 2 O 2 , (b) O 3 , CO, (c) NO and PAN obtained in KPP and KINAL.
Author Contributions: L.C. initiated the study presented in this paper, and instructed all the simulation scenarios. S.L., Z.Y. and L.C. performed the computations, plotted the results, and wrote the manuscript together. M.G. was involved in the discussion of the computational results and revised the manuscript. All the authors listed have read the final version of the manuscript and approved the submission.