Investigation of the wet removal rate of black carbon in East Asia: validation of a below- and in-cloud wet removal scheme in FLEXPART v10.4

Abstract. Understanding the global distribution of atmospheric black carbon (BC) is essential to unveil its climatic effect. However, there are still large uncertainties regarding the simulation of BC transport due to inadequate information about the removal process. We accessed the wet removal rate of BC in East Asia based on long-term measurements over the 2010–2016 period at three representative background sites (Baengnyeong and Gosan in South Korea and Noto in Japan). The average wet removal rate, represented by transport efficiency (TE), i.e. the fraction of undeposited BC particles during transport, was estimated as 0.73 in East Asia from 2010 to 2016. According to accumulated precipitation along trajectory, TE was lower in East and North China, where the industrial sector (thin-coated) is dominant; in contrast, that in South Korea and Japan showed higher values due to the transport sector (thick-coated), with emissions mainly from diesel vehicles. By the same token, TE in winter and summer showed the highest and lowest values, respectively, depending on the dominant emission sectors, such as house heating (thick-coated) and industry. The average half-life and e-folding lifetime of BC were 2.8 and 7.1 days, respectively, similar to previous studies, but those values differed according to the geographical location and meteorological conditions of each site. Next, by comparing TE from the FLEXible PARTicle (FLEXPART) Lagrangian transport model (version 10.4), we diagnosed the scavenging coefficients (s−1) of the below- and in-cloud scavenging scheme implemented in FLEXPART. The overall median TE from FLEXPART (0.91) was overestimated compared to the measured value, implying underestimation of wet scavenging coefficients in the model simulation. The median of the below-cloud scavenging coefficient showed a lower value than that calculated from FLEXPART, by a factor of 1.7. On the other hand, the overall median of the FLEXPART in-cloud scavenging coefficients was highly underestimated by 1 order of magnitude compared to the measured value. From the analysis of artificial neural networks, the convective available potential energy, which is well known as an indicator of vertical instability, should be considered in the in-cloud scavenging process to improve the representative regional difference in BC wet scavenging over East Asia. For the first time, this study suggested an effective and straightforward evaluation method for wet scavenging schemes (both below- and in-cloud) by introducing TE along with excluding effects from the inaccurate emission inventories.



Introduction 40
Black carbon (BC) is the most significant light-absorbing aerosol that can cause positive radiative forcing on climate change 41 (Winiger et al., 2016;Myhre et al., 2013;Bond et al., 2013;Emerson et al., 2018). However, state-of-the-art models still have 42 a limitation in evaluating the direct radiative forcing of BC because of the large model uncertainties in simulating BC 43 concentrations (Xu et al., 2019;Bond et al., 2013;Samset et al., 2014;Wang et al., 2014a). This can partly be attributed to the 44 following three reasons: inaccurate bottom-up emission inventory, the complexity of BC hygroscopicity, and an imprecise 45 dry/wet deposition scheme. First, because BC is mainly contributed by scattered emission sources, the uncertainty of BC 46 emission rates is large compared to other species whose emissions are dominated by large sources (Zheng et al., 2018). Second, 47 Figure 1c reveals the geographical distribution for the mean BC mass of identified potential emission regions, indicating that 135 this approach was appropriate because of good spatial coverage over East Asia, including East China, a major emission source 136 for BC. 137

Transport efficiency (TE) 138
The TE of BC is defined as the ratio of the BC and CO concentrations measured at the receptor site to that anticipated if 139 there was no wet removal during transport (i.e., APT is zero). Thus, the TE of the airmass was calculated by eq. (1), 140 where delta (Δ) indicates the difference between BC and CO concentrations and their baseline concentrations (Moteki et al., 142 2012;Oshima et al., 2012;Kanaya et al., 2016). The baseline CO was estimated as a 14-day moving 5th percentile from the 143 observed CO mixing ratio, but the BC baseline was regarded as zero because the atmospheric lifetime of BC is known as 144 several days, which is much shorter than that of CO (1−2 months). [ΔBC/ΔCO]APT=0 indicated the regional median value of 145 ΔBC/ΔCO under dry conditions implying the original emission ratio. In our previous work, we successfully elucidated that 146 [ΔBC/ΔCO]APT=0 depends on the regional characteristics of the energy consumption types (Kanaya et al., 2016;Choi et al., 147 2020). The decrease of the ratio with APT, [ΔBC/ΔCO]APT>0, was related to BC-specific removal due to wet scavenging 148 processes and thus the TE was effective indicator to investigate the wet removal process. Although TE is also affected by dry 149 https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. 5 deposition, but the effect of dry deposition could be negligible because dry deposition velocities were much lower than the 150 default setting in the global models (Choi et al., 2020). 151

FLEXPART model 152
To compare the wet removal rates between the model simulation and measured values, the FLEXPART v10.4 was used to 153 simulate BC wet scavenging over East Asia using the backward mode. Detailed information for the FLEXPART is readily 154 found in the literature (e.g., Stohl et al., 2005); thus, we only briefly describe the information here. The FLEXPART version 155 10.4 was the official version to allow turning on the wet scavenging module in the backward simulation mode 156 (https://www.flexpart.eu/downloads, obtained 10 October 2019). The equation and detailed description for the below-and in-157 cloud scavenging scheme are explained in Pisso et al. (2019) and Grythe et al. (2017). The FLEXPART model was executed 158 with operational reanalysis meteorological data from the ECMWF ERA-Interim at a spatial resolution of 1°×1° with 60 full 159 vertical levels. Temporally, ECMWF has a resolution of 3 h, with 6 h analysis and 3 h forecast time steps. The period and 160 daily frequency of simulation were the same as those of the HYSPLIT model (past 72 h and four times, respectively). It should 161 be noted that chemistry and microphysics could not be resolved by the FLEXPART. The FLEXPART model, therefore, ignores 162 the aging process (from hydrophobic to hydrophilic state changes and size changes of BC) and assumes that all BC particles 163 are aged hydrophilic particles. The logarithmic size distribution of BC with a mean diameter of 0.16 μm and a standard 164 deviation of 1.84, in accordance with measurement in Japan, was used (Miyakawa et al., 2017). A total of 10 4 particles were 165 randomly released at 500 m from each receptor site during 1 h when the measurement data were existed. To validate the wet 166 scavenging scheme in FLEXPART by comparison with the measured TE value, the wet scavenging coefficients for below-and 167 in-clouds were extracted from FLEXPART to calculate TE. 168 3 Results 169

Overall variation of transport efficiency (TE) 170
Figure 2 shows that measured [ΔBC/ΔCO]APT=0 (left panel) and TE variations (right panel) depend on APT and the 171 measurement sites. The overall median [ΔBC/ΔCO]APT=0 was 6.4 ng m −3 ppb −1 , which converged from Baengnyeong (6.2 ng 172 m −3 ppb −1 ), Gosan (6.5 ng m −3 ppb −1 ) and Noto (6.7 ng m −3 ppb −1 ), indicating that TE is characterized with a regional 173 [ΔBC/ΔCO]APT=0 per site. We divided APT into 9 range bins and applied exponential fitting equations to quantify the wet 174 removal process. Among NAPT > 0 (total number of data points when APT > 0 mm), only the data point fraction in each bin to 175 NAPT > 0 ≥ 2% was considered to secure the statistic. It should be noted that we found the relationship between TE and APT by 176 using the stretched exponential decay (SED) equation,

instead of the widely used equation, A-177
B×log(APT), because the coefficients of determination (R 2 ) was improved up to 0.981 though TE values from three sites were 178 used (Table 1). This fitting equation is normally used to describe below-cloud scavenging, whereas wet removal of BC is 179 generally believed to be dominated by in-cloud rather than below-cloud processes because of the small size of BC-containing 180 particles. Therefore, the equations should contain both below-and in-cloud scavenging effects. The parameters A1 (0.269 ± 181 0.039) and A2 (0.385 ± 0.035) of the overall fitting were higher and lower, respectively, than the derived equation from the 182 Fukue site (A1 = 0.109 and A2 = 0.68) (Kanaya et al., 2016). It can be easily deduced that the wet removal effect at the three 183 sites was initially more effective than that at Fukue, but the wet removal effect at Fukue gradually accelerated as the APT 184 increased. In particular, the A2 value is important for calculating the TE of BC for long-range transport, e.g., toward the Arctic 185 https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. 6 (Kanaya et al., 2016;Zhu et al., 2019), because A2 determines the magnitude of the wet removal efficiency according to APT. 186 Thus, the newly obtained SED equation indicates that more BC will be transported to the Arctic region than previously reported. 187 The decreasing pattern of median TE for Baengnyeong did not closely follow the overall SED and had a much lower R 2 188 (0.77), indicating that the wet removal process at Baengnyeong could not simply be expressed by APT. In contrast, the R 2 of 189 Gosan and Noto were sufficiently high to represent the wet removal characteristics. The aging process due to different traveling 190 times might be one of the reasons. Because long-range transported BC has a larger core diameter than BC from local sources 191 (Lamb et al., 2018;Ueda et al., 2016), these larger BC cores are preferentially removed via the wet scavenging process (Moteki 192 et al., 2012). However, previous studies reported that the mass median diameter (MMD) of BC at Baengnyeong, Gosan,and 193 Noto in spring were 218, 196, and 200 nm (Oh et al., 2015;Ueda et al., 2016;Oh et al., 2014), respectively, indicating much 194 more aging compared with local emissions in Seoul, South Korea (180 nm) and Tokyo, Japan (163 nm) (Park et al., 2019;195 Ohata et al., 2019). Moreover, there were no significant differences in the mean traveling times for the airmass (when APT > 196 0) arriving at the three sites (37.9, 39.0, and 37.8 h for Baengnyeong, Gosan, and Noto, respectively), indicating that the 197 difference in the level of the BC aging process might be negligible. 198 The difference in the wet removal rate among measurement sites could be partly explained by the difference in meteorology. 199 The monthly mean meteorological parameters indicated that Baengnyeong has characteristics of low precipitation (80.6 mm), 200 cloud cover (0.57), total column cloud water (0.06 kg m −2 ), and high cloud bottom height (2.5 km) compared to other sites, 201 suggesting the lower exposure time to both below-and/or in-cloud condition during the transportation (Figure 3). In contrast, 202 the SED fittings for both Gosan and Noto showed similar ranges of high precipitation (127 and 174 mm), total cloud cover 203 (0.65 and 0.64), and total column cloud water (0.09 and 0.12 kg m −2 ) but low cloud bottom height (1.9 and 2.0 km), respectively. 204 In addition, the different BC coating thicknesses according to the emission source and fuel types could also contribute to the 205 site difference of the wet removal rate, which will be further discussed in section 3.2. 206 Using the overall SED fitting equation, TE at 0.5 (TE=0.5) and e-folding (TE=1/e) could be reached when the APT values 207 were 11.7 and 30.2 mm, respectively (Table 1). Similar to the SED results, Baengnyeong needed much higher precipitation of 208 70.9 and 202 mm to reach TE=0.5 and TE=1/e, respectively, but the other sites showed lower APTs of 16.4 mm and 42.3 mm 209 for Gosan and 8.0 mm and 20.3 mm for Noto, respectively. Considering the annual mean precipitation at the three sites (1542 210 mm), it took 2.8 and 7.1 days to reach TE=0.5 and TE=1/e, respectively. Kanaya et al. (2016) reported a similar half-life and 211 shorter e-folding lifetime for BC at Fukue (2.3 ± 1.0 and 4.0 ± 1.0 days, respectively), calculated from the 15.0 ± 3.2 mm and 212 25.5± 6.1 mm of APT to reach TE=0.5 and TE=1/e, respectively, along with annual precipitation (2335 mm). This calculated 213 e-folding lifetime in East Asia was much shorter than 16.0 days for the global model (Grythe et al., 2017). 214 Based on a similar approach over the Yellow Sea using an aircraft-borne single particle soot photometer (SP2) during the 215 A-FORCE campaign , attaining TE=0.5 required different magnitudes of APT depending on not only the 216 airmass origin but also the altitude. These authors also reported that the TE of northern China was higher than that of southern 217 China regardless of altitude. Therefore, in the next section, we will further investigate why the difference in halving or e-218 folding lifetimes depends on region and season by analyzing the difference in the origin of airmasses and the seasonal variation 219 of BC emission sources.  (Table 1). A similar tendency of R 2 , TE=0.5 also showed different APTs, i.e., higher in East and North China and 224 lower in other regions. The regional difference in wet removal efficiency can partly be attributed to the following reasons. 225 First, the transport pathway of airmasses from East and North China could be less exposed to in-cloud scavenging than other 226 regions because the most of potential emission source in East and North China is located over 30°N (Figure 1c), which has 227 low cloud cover and water contents along with high cloud bottom heights ( Figure 3). Although the amount of APT was similar 228 to other regions, it was mostly composed of below-cloud scavenging, therefore, the wet removal efficiency should be lower 229 than the dominant in-cloud scavenging region. Second, the difference in the coating thickness of BC particles, depending on 230 the emission sectors, could be a major factor causing the difference in the wet removal efficiency because thickly-coated BC 231 particles are much easier to remove by wet scavenging than less coated and/or freshly emitted BC (Miyakawa et al., 2017). 232 Typically, BC emitted from industrial regions, transport from diesel vehicles, and domestic sectors has characteristics of weakly, 233 moderate, and strongly coated BC, respectively (Han et al., 2019;Liu et al., 2019), based on insignificant differences in the 234 MMD of BC from those emission sectors (190 -200 nm). This result coincided with the major emission sector of the REAS 235 emission inventory in East and North China and North Korea (~57.5% emitted from industrial sectors) compared to other sites 236 (12% − 39%). In contrast, Northeast China showed low APT for reaching TE=0.5 and TE=1/e because the dominant BC 237 emission sector was residential sector (48.3%) which has a thickly coated characteristic. BC from South Korea and Japan 238 reached TE=0.5 and TE=1/e with a small amount of APT because moderately coated BC was mostly emitted from the transport 239 sector (73.4%), mainly from diesel vehicles. It should be noted that the dominant emission sectors of industry (for East and 240 North China and North Korea) or transport sectors (South Korea and Japan) were also confirmed by the Emission Database 241 for Global Atmospheric Research (EDGAR) in 2010 and MIX in 2010 (Li et al., 2017;Crippa et al., 2018). 242 In case of seasonal variation in TE, the decreasing magnitude of TE was obviously emphasized in fall and winter, which 243 was much steeper than that in spring and summer ( Figure 4b). This tendency was reflected in the effect of the residential sector, 244 which has thickly coated BC, which increased due to house heating as the temperature decreased. In contrast to winter, the 245 APT for reaching TE=0.5 in spring and summer was the highest among the seasons. This might be caused by the increasing 246 fraction of BC from the industrial sector in China while decreasing emissions from residential sectors (Kurokawa et al., 2013). 247

Comparison of measured and FLEXPART-simulated TE 248
In this section, by extracting the wet scavenging coefficients (Λ; s -1 ) from the FLEXPART simulation, the difference in TE 249 between the measured and simulated values was investigated. The scavenging coefficient (Λ; s -1 ) is defined as the rate of 250 aerosol washout and/or rainout due to the wet removal process. The TE value based on measurements and FLEXPART can be 251 expressed by multiplying each TE (1 -removal rate) of serial grid cells as in eq. (2), 252 where ηn indicates the removal rate in the nth grid cell and is expressed as eq. (3), 254 where t and fg indicate the residence time and fraction for the subgrid in a grid cell, respectively. Because the precipitation is 256 not uniform in a single grid cell, fg accounts for the variability of precipitation in a grid cell in FLEXPART. fg is a function of 257 https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. 8 large-scale and convective precipitation, as described in Stohl et al. (2005). Although the grid resolution of the input 258 meteorological data for the HYSPLIT model (0.25°×0.25°) is much finer than that for FLEXPART (1°×1°), we assumed the 259 same potential emission region as the HYSPLIT model for calculating TE because there was no significant difference in the 260 airmass pathway between the two outputs. 261 The overall median value of measured TE was 0.72, and Baengnyeong showed the highest (0.88), followed by Gosan (0.70) 262 and Noto (0.68) due to reasons explained in the previous sections. In comparison, the overall median value of FLEXPART TE 263 (0.91) was much higher than the measured TE, indicating that the wet scavenging coefficients in the FLEXPART scheme were 264 significantly underestimated. Moreover, the difference in FLEXPART TE depending on the measurement sites (0.95 for 265 Baengnyeong, 0.94 for Gosan, and 0.87 for Noto) was not large as the measured TE, suggesting that the regional difference in 266 meteorological variables was relatively normalized and that the influence of other variables, which were not considered in the 267 wet scavenging scheme, might be excluded in the calculation. Meanwhile, it is difficult to capture the local variation from 268 coarse grid sizes, despite the airmass transport pathway between the two models being similar, because the key variables for 269 determining the wet scavenging coefficient (such as precipitation and cloud cover) could have a large local variability. In 270 addition, this approach still had a limitation in determining whether the overestimation of TE was resulting from the below-or 271 in-cloud scavenging processes. Nevertheless, with similar rationale, further comparison of measured and simulated scavenging 272 processes would provide information to better represent wet removal schemes. 273

Below-cloud scavenging efficiency (Λbelow) 274
From this section, we aimed to investigate the below-and in-cloud scavenging in detail by discriminating the representative 275 cases according to cloud information from the ERA5 pressure level data to overcome the limitation of the local variability of 276 meteorological input variables. By considering the vertical height of the airmass from the HYSPLIT model and cloud 277 information from ERA5, we successfully distinguished the dominant cases for below-cloud (no residence time within the cloud) 278 and in-cloud (no residence time below the cloud) cases when precipitation ≥ 0.01 mm hr −1 . The median TE and residence time 279 for only in-cloud cases (0.72 and ~7,200 h) were much lower and longer, respectively, than those for only below-cloud cases 280 (0.89 and ~5,100 h), indicating that most BC particles were effectively removed via the in-cloud scavenging process (Table 2). 281 In the case of below-cloud scavenging, the deviation of TE from unity could be simply converted to the scavenging coefficient 282 (Λbelow) by considering the precipitation intensity, raindrop size, aerosol size, and residence time in a grid cell. Because many 283 studies have made an effort to parameterize Λbelow using observation data and/or the theoretical calculation (Xu et al., 2017;284 Wang et al., 2014b;Feng, 2007), we also parameterized this coefficient using a simplified method by following the scheme of 285 below-cloud scavenging in FLEXPART v10.4 (Laakso et al., 2003), which only considers the precipitation rate and aerosol 286 size. Assuming a BC size ~200 nm, TE for below-cloud can be expressed using equations (2) and (3)

by substituting Λ with 287
Λbelow, which depends only on the precipitation rate in the subgrid cell (Itotal; the ratio of precipitation to fg). Because Λbelow can 288 be determined by constraining the proportion to the summation of Itotal, hourly Λbelow from the sequential grid cell in a single 289 case can easily be obtained by minimizing χ 2 , (TEmeasured − TEcalculated) 2 when χ 2 < 0.1. This was conducted using an R function, 290 optimization (optim; https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html), included in the standard R package 291 "stats". 292 Figure 5a indicates the empirical cumulative density function for the measured Λbelow from 869 cases. Although a substantial 293 fraction of Λbelow was close to zero (or negative), the median Λbelow was significantly different from zero and also positive 294 (7.9×10 -6 s -1 ), with an interquartile range of -1.7×10 -5 s -1 to 5.3×10 -5 s -1 . Negative Λbelow values have been reported in previous 295 https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. 9 studies (Laakso et al., 2003;Pryor et al., 2016;Zikova and Zdimal, 2016); therefore, we assumed that these negative values 296 reflected the uncertainty in measurements and/or inclusion of BC, which might be continuously supplemented in airmasses. 297 As the threshold of Itotal increased from 0.01 (all cases) to 0.2 mm hr -1 (median), Λbelow values were increased by a factor of 2.5 298 to 2.0×10 -5 s -1 (-2.5 ×10 -5 s -1 to 9.0 ×10 -5 s -1 ). Using these obvious increasing tendencies of Λbelow according to Itotal, we 299 determined the empirical fitting equation by investigating the relationship between median Λbelow and each bin of Itotal. Figure  300 5b indicates Λbelow as a function of Itotal by allocation to 11 logarithmic bins. As the estimated Itotal bins covered the Itotal ranges, 301 0.03 to 2.0 mm hr -1 (5 th percentile to 95 th percentile), this exponential fitting equation (A × Itotal B ) could be representative for 302 below-cloud scavenging over East Asia. The constant A and exponent B with a 95% confidence interval were 2.0 × 10 -5 (1.9 303 -2.2 × 10 -5 ) and 0.54 (0.46 -0.64), respectively. Instead of the SED equation shown in Figure 2, we chose the exponential 304 fitting equation because of its higher R 2 (0.973) compared to that from SED fitting (0.903), as well as being widely used in 305 previous studies. 306 constant with increasing Itotal because the mean absolute MFBs were slightly increased from 1.4 to 1.6. It should be noted that 311 Λbelow from Laakso et al. (2003), which is the default scheme for below-cloud scavenging in the FLEXPART model version 312 10 or higher (Grythe et al., 2017), showed fairly good agreement with our Λbelow among the reported values (mean absolute 313 MFBs was 0.68). MFB was positive at low Itotal, but the tendency was opposite based on Itotal ~ 0.1 mm hr -1 , suggesting that 314 Λbelow might be converged within a similar range when we consider the range of Itotal. Although Λbelow from Laakso et al. (2003) 315 showed good agreement with our results, the median Λbelow (6.6 × 10 -6 s -1 ) was overestimated compared to our estimation (4.0 316 × 10 -6 s -1 ), by a factor of 1.7 when we recalculated the only below-cloud cases. The MFBs from other schemes were too high 317 or low to declare reasonable results. For example, the Λbelow of secondary ions in Beijing (Xu et al., 2017) had the highest MFB 318 (1.68), and although the diameter ranges were larger (~ 500 nm) than those of BC, the effect of differences in diameter might 319 be negligible. 320

In-cloud scavenging coefficient (Λin) 321
Compared to Λbelow, the calculation of Λin is much more complicated because many factors can influence the in-cloud 322 scavenging process, such as precipitation, total cloud cover (TCC), the specific cloud total water content (CTWC) and so on. 323 The detailed description for the complicated equation for Λin in FLEXPART v10 is presented in Grythe et al. (2017) where icr and Fnuc are the cloud water replenishment factor (6.2; default value) and the nucleation efficiency, respectively. It 327 should be mentioned that Λin was calculated by following the FLEXPART scheme using the ERA5 meteorological data 328 (0.25º×0.25º) instead of the FLEXPART simulation (1º×1º) to match the grid size of the input data with the HYSPLIT backward 329 trajectory. Among the 769 cases for in-cloud cases, equations (2) and (3) were also used to calculate TE for only in-cloud cases 330 by substituting Λ with calculated Λin. Unlike the hourly measured Λbelow calculated by optimization, the only overall median 331 https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. Λin (Λin*) for in-cloud cases was calculated using equation (3) because Λin cannot be constrained by a specific variable. 332 The FLEXPART Λin* (7.28 × 10 -6 s -1 ) was underestimated by 1 order of magnitude compared to our estimated Λin* (8.06 × 333 10 -5 s -1 ). When TE from FLEXPART for in-cloud cases (all cases) was recalculated by considering a ten (five) times higher 334 Λin, the median TE was 0.73 (0.79), which was much close to the measured TE (0.72). Although the grid size of input 335 meteorological data for two approaches did not match, the underestimation of the in-cloud scavenging scheme in FLEXPART 336 was confirmed. Grythe et al. (2017) reported an overestimation of observed BC (a factor of 1.68) due to the inaccurate emission 337 source rather than the underestimated in-cloud removal efficiencies. Although the effect of BC particle dispersion to adjacent 338 grid cells was neglected in our approach, the underestimation of in-cloud scavenging coefficients was obvious because the 339 accuracy of the emission inventory did not affect the estimated Λin*. Looking more closely into the sites, the FLEXPART Λin* 340 at Noto was remarkably underestimated by 1 order of magnitude, followed by Gosan (~90%) and Baengnyeong (~43%), 341 similar to the order of the wet removal efficiency. It should be noted that the coefficient of variation (CV; standard deviation 342 divided by the mean) of FLEXPART Λin* was much lower (0.23) than the measured Λin* (0.78), indicating that FLEXPART 343 Λin* did not accurately represent the actual regional difference in the real world. Among the input meteorological variables in 344 equation (4), the CV of Itotal was the highest as 0.22, which was similar to the CV of FLEXPART Λin*, followed by CTWC 345 (0.08), fg (0.03), and TCC (0.02), suggesting that the difference in FLEXPART Λin* could be partially explained by Itotal rather 346 than other variables. Among the meteorological variables that were not considered in equation (4), the convective available 347 potential energy (CAPE), which is well known as an indicator of vertical instability (Mori et al., 2014), had the highest CV of 348 0.31. 349 We employed an artificial neural network (ANN) method to compare the importance of CAPE with other considered input 350 meteorological variables for determining the hourly Λin, not Λin*. We applied a stricter selection for in-cloud cases that only 351 when in-cloud scavenging occurred less than three times (i.e., three cells) in a single case, regardless of the number of below-352 cloud occurrences. Because the effect of below-cloud scavenging was successfully excluded from the TE using the derived 353 equation for Λbelow in the previous section, the Λin in less than three in-cloud cases can also be calculated by optimization based 354 on the remaining TE. We applied a threshold of three cases here because the number of data (230 cases) was sufficient to 355 conduct statistical analysis, while the optimization uncertainty could be reduced to its minimum. The ANN model was trained 356 using the six meteorological variables (CAPE, CTWC, fg, Fnuc, Itotal, and TCC), and all variables were normalized by the 357 minimum and maximum of each variable ([x-min(x)]/[max(x)-min(x)]). To determine optimal node numbers in the hidden 358 layer, we applied a 'caret' package of the R function that contain several sets of machine learning modes and validation tools 359 (https://cran.r-project.org/web/packages/caret/caret.pdf) and adopted a method from the 'neuralnet' package that is fit for a 360 multi-hidden layer. By varying the 'size' (node number) from 5 to 20 and using k-fold cross validation, the selected cases were 361 randomly divided 3:1 into training (172 data points) and validation data (58 data points). Garson's algorithm in the 362 "NeuralNetTools" package was used to identify the relative importance of six input variables in the final neural network 363 (Garson, 1991). The model's performance was assessed in these independent validation data by calculating the root mean 364 squared error. The optimal number of nodes in the hidden layer was 12 (Figure 7a). 365 Figure 7b shows the relative importance of input variables for calculating Λin using Garson's algorithm. The most important 366 input variable was CAPE, with a value of 35%, followed by CTWC, Itotal, and so on, confirming that CAPE should be 367 considered in the Λin calculation. Typically, enhancing wet removal by convective clouds successfully reduces the aloft BC 368 concentration in the free troposphere (Koch et al., 2009). Therefore, convective process is important in tropical regions but has 369 a slightly lower impact at mid-latitudes (Luo et al., 2019;Grythe et al., 2017;Xu et al., 2019). Moreover, previous studies 370 pointed out convective scavenging to be a key parameter in determining the BC concentration in model simulation (Lund et 371 al., 2017;Xu et al., 2019) and the role of wet removal by convective clouds might be significant when most airmasses travel 372 above the planetary boundary layer. Unfortunately, the current version of FLEXPART does not implement convective 373 scavenging (Philipp and Seibert, 2018), which could be a plausible reason for the underestimation of FLEXPART Λin. Although 374 the relative importance of each variable cannot be parameterized to calculate Λin, this approach highlights that CAPE is one 375 of the key factors for determining Λin over East Asia. In the future, more information might be required to evaluate the in-376 cloud scavenging scheme using Weather Research and Forecasting (WRF)-FLEXPART at a higher resolution in further studies, 377 since a 0.25º grid size is still not sufficient to reproduce convective clouds (typically 10 km or less). 378

Conclusions 379
The wet removal rates and scavenging coefficients for BC were investigated by the term of ΔBC/ΔCO ratios from long-380 term, best-effort observations at three remote sites in East Asia (Baengnyeong and Gosan in South Korea and Noto in Japan). 381 Combined with backward trajectories covering the past 72 h, accumulated precipitation along trajectories (APT), and transport 382 efficiency (TE; [ΔBC/ΔCO]APT>0/[ΔBC/ΔCO]APT=0), the assessment of BC wet removal efficiency was conducted as an aspect 383 of the pathway of trajectories, including the successful identification of below-and in-cloud cases. The overall wet removal 384 rates as a function of APT, the half-life and e-folding lifetime were similar to those of previous studies but showed large 385 regional differences depending on the measurement sites. The difference in the wet removal rate, depending on the 386 measurement site, can be explained by the different meteorological conditions, such as the precipitation rate, cloud cover, total 387 column cloud water, and cloud bottom height. Moreover, the difference in regional or seasonal wet removal rates could be 388 explained by the different coating thicknesses according to the BC emission sources (thin-and thick-coated BC from the 389 industrial and residential sectors, respectively) because the thick-coated BC particles are preferentially removed due to cloud 390 processes. By discriminating below-and in-cloud dominant cases according to cloud vertical information from ERA5 pressure 391 level data, scavenging coefficients for below-cloud (Λbelow) and in-cloud (Λin * ) were simply converted from the measured TE 392 values. The Λbelow from the FLEXPART scheme was overestimated by a factor of 1.7 compared to the measured Λbelow, although 393 the measured Λbelow showed good agreement with the below-cloud scheme in FLEXPART among the reported scavenging 394 coefficients. In contrast to Λbelow, FLEXPART Λin * was highly underestimated by 1 order of magnitude compared to measured 395 Λin * , suggesting that the current in-cloud scavenging scheme did not represent regional variability. By diagnosing the relative 396 importance of the input variables using the artificial neuron network (ANN) method, we found that the convective available 397 potential energy (CAPE), which is an indicator of vertical instability, should be considered to improve the in-cloud scavenging 398 scheme because convective scavenging could be regarded as a key parameter for determining the accurate BC concentration 399 in a model. This study could contribute not only to improving the below-cloud scavenging scheme implemented in a model, 400 especially FLEXPART, but also to providing evidence for complementary in-cloud scavenging schemes by considering the 401 convective scavenging process. For the first time, these results suggest a novel and straightforward approach to evaluating the 402 wet scavenging scheme in various models and to enhancing the understanding of BC behavior by excluding the effect of 403 inaccurate emission inventory. 404 https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. Figure 1. (a) The location of three measurement sites (Baengnyeong, Gosan, and Noto) and the black carbon (BC) emission rate (ton year -1 ) over East Asia from the Regional Emission inventory in ASia (REAS) version 2.1 (Kurokawa et al., 2013). (b) Illustration of residence time calculated based on the HYSPLIT backward trajectory that passed over a single grid cell (see details in the manuscript). (c) The spatial distribution of the mean BC mass in the potential emission region, which is the highest BC mass grid of each trajectory. The BC mass was obtained by multiplying (a) the emission rates and (b) the residence time.
https://doi.org/10.5194/acp-2020-402 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License.     The relative importance of six input meteorological variables used for calculating in-cloud scavenging coefficients in the FLEXPART model (except for CAPE) using Garson's algorithm implemented in the 'NeuralNetTools' package in R. CAPE, CTWC, Itotal, fg, TCC, and Fnuc represent the convective available potential energy, specific cloud total water content, precipitation rate, fraction of a subgrid in a grid cell (see manuscript for details), total cloud cover, and nucleation efficiency, respectively.