ChAP 1.0: A stationary tropospheric sulphur cycle for Earth system models of intermediate complexity

. A stationary, computationally efﬁcient scheme ChAP-1.0 (Chemical and Aerosol Processes, version 1.0) for the sulphur cycle in the troposphere is developed. This scheme is designed for Earth system models of intermediate complexity (EMICs). The scheme accounts for sulphur dioxide emissions into the atmosphere, its deposition to the surface, oxidation to sulphates, and dry and wet deposition of sulphates on the surface. The calculations with the scheme are performed forced by anthropogenic emissions of sulphur dioxide into the atmosphere for 1850-2000 adopted from the CMIP5 dataset and by 5 the ERA-Interim meteorology assuming that natural sources of sulphur into the atmosphere remain unchanged during this period. The ChAP output is compared to changes of the tropospheric sulphur cycle simulations: with the CMIP5 data, with the IPCC TAR ensemble, and with the ACCMIP phase II simulations. In addition, in regions of strong anthropogenic sulphur pollution, ChAP results are compared to other data, such as the CAMS reanalysis, EMEP MSC-W, and with individual model simulations. Our model reasonably reproduces characteristics of the tropospheric sulphur cycle known from these information 10 sources. In our scheme, about half of the emitted sulphur dioxide is deposited to the surface and the rest in oxidised into sulphates. In turn, sulphates are mostly removed from the atmosphere by wet deposition. The lifetime of the sulphur dioxide and sulphates in the atmosphere is close to 1 day and 5 days, respectively. The limitation

The latter model also implements a very simple atmospheric sulphur cycle scheme, in which sulphate burden per unit area is related to their precursor emissions at the same grid cell by using a prescribed coefficient, which, in turn, is related to atmospheric lifetimes of sulphates and their precursors taking into account that part of the emitted precursors are deposited before they are oxidised into sulphates. No horizontal transport of sulphates and their chemical precursors are :: is allowed for.
This approach is reasonable for Climber-2 with its very coarse horizontal resolution (10 • by latitude and 51.3 • on longitude, compounds) with the vertical scales H Y which is of the order of 1 km (Jaenicke, 1993;Warneck, 2000). The latter leads to the 90 relation between the near-surface mass concentrations q Y,s and the total burden B Y of substance Y per unit area: This allows us to integrate (1) over vertical coordinate, formally from the surface up to the infinity. To simplify the setup, we neglect the dependence of horizontal velocity on the vertical coordinate. The resulting equations read 95 where Y ∈ {SO 2 , SO 4 }, (3) is similar to (1) with two important differences: now U is two-dimensional (we choose it at a representative altitude), and A Y represents only horizontal diffusion.
Further, we assume that major chemical reactions follow the common first-order kinetics relative to the source compounds. 100 In a similar fashion, we assume that sink terms are proportional to the respective burdens. All this leads to Here k's stand either for the respective reaction rate constants or for the loss rate coefficients.
Based on (Warneck, 2000) and on individual model simulations summarised in (Houghton et al., 2001, their Table 5.5), we neglect the following terms in (1)

105
both non-stationary terms ∂B Y /∂t, chemical SO 2 production in the atmosphere: R SO2,prod = 0 (this assumption basically removes part of the natural sources of sulphur dioxide, e.g., the DMS oxidation), gas-phase sulphur dioxide oxidation: R gas = 0, wet deposition of sulphur dioxide: D SO2,wet = 0 (or, equivalently, k SO2,wet = 0), 110 This and the previous assumption sets make Eqs. (3) linear with respect to prognostic variables. For time being, we additionally drop diffusion terms A Y . Thus, Eqs. (3) are reduced to Here, k SO2 = k in−cl + k SO2,dry , k SO4 = k SO4,dry + k SO4,wet , and E SO2 is SO 2 emission rate per unit area.
To get a guide :: for ::: the :::::: further, consider an one-dimensional problem with U = |u| = const and the emission localised in the 115 interval 0 ≤ x ≤ L, where x is coordinate in the direction of U , and L is the horizontal source size (Jacob, 2000). Its solution are shown in supplementary Fig. S1. The averaged over the 0 ≤ x ≤ L domain the solution of the first equation (5) reads where γ Y = k Y L/U . Downwind of the emission region this solution reads 120 Eq. (7) allows, in particular, to estimate the horizontal scale of influence for this grid cell source as where T Y = k −1 Y is the lifetime of species Y in the atmosphere. An estimate of T SO2 may be obtained from Eq. (8) by using the ERA-Interim data (Dee et al., 2011). For this dataset, the typical values of zonal u and meridional v velocities in the lower troposphere in the middle latitudes, where the SO 2 pollution 125 is most marked, are up to 7 m s −1 and up to 5 m s −1 (Fig. S2). Because the typical SO 2 lifetime is around 1-2 days (Warneck, 2000;Surkova, 2002;Houghton et al., 2001, their Table 5.5), we can estimate L SO2 ∼ 10 3 km (somewhat larger in the zonal direction and somewhat smaller in the meridional one). This corresponds to few grid cells provided that the grid cell size is several hundred kilometres.
Typical horizontal scale for SO 4 , L SO4 , may be estimated in a similar way. Assuming that the transport velocity is of the 130 same order of magnitude as it was used for SO 2 advection (this assumption is justified by similar depths of the atmospheric layers -around 1.5−2 km (Warneck, 2000) and taking 4−6 days as a typical value for T SO4 (Warneck, 2000;Surkova, 2002;Houghton et al., 2001, their Table 5.5), we estimate L SO4 ∼ 5 × L SO2 .

Horizontal transport
Taking into account just obtained estimates for L SO2 and L SO4 , we may construct the transport scheme as follows. 135 At first, we solve the first Eq. (5). Because it is linear with respect to B SO2 , we may consider the model grid as an array of non-interacting sulphur sources numbered as j = 1, 2, . . .. To reduce the computational burden, we consider only grid cells with E (j) SO2 ≥ E SO2,min and set E SO2,min to sufficiently small, empirically chosen value (Table 1). At the source grid cell, the burden B (j) SO2 (ρ j ) is calculated by using Eq. (6) with Y = SO 2 and with γ SO2 = k SO2 (∆x/u + ∆y/v). Here ρ j is the horizontal coordinates of the grid cell corresponding to the source j, E (j) SO2 = E SO2 (ρ j ), ∆y = R E ∆φ, ∆x = R E cos φ∆λ, 140 R E is Earth radius, φ is latitude, ∆φ and ∆λ are grid cell sizes in latitudinal and longitudinal directions, respectively.
The difference between E (j) SO2 and k SO2 (ρ j )B (j) SO2 (ρ j ) is transported out of the source cell by advection: This flux is partitioned into zonal and meridional components in proportion to the corresponding wind component and to the geometric size of the corresponding boundary of the cell: The direction of each F is determined by the direction of zonal and meridional wind respectively. Then we loop in the zonal direction and calculate SO 2 burdens, related to the source j, as well as corresponding chemical and depositional losses, and fluxes out of the cell by using Eqs. (6), (9), and (10) (Fig. 1). In each cell i included in this loop, zonal flux from the previous zonal cell F loop is stopped in the cell with the number I (j) , at which any of the conditions is met: either zonal wind u(ρ I (j) ) changes sign relative to u(ρ j ) or the whole latitudinal circle is looped over.
In the stopping cell, F The SO 2 mass, which is advected from the source cell and from the each of the looped-over cells, is transported to the 155 respective neighbour cell either to the north or to the south depending on sign of v(ρ i ). No advection of sulphur dioxide mass from this meridional neighbour cell is allowed, and the SO 2 burden is calculated assuming the balance between meridional mass inflow and chemical loss ( Fig. 1).
At the next step, the SO 2 burden in each grid cell is obtained by summing over all grid-cell sources: where SM O is the just described smoothing operator. We set n smo equal to 5 in the contemporary implementation.
It is easy to show that, by construction, for each grid-cell source E (j) where i stands for the set of cells over which the zonal loop is performed together with their meridional neighbours (or, in other words, the set of all coloured cells in Fig. 1). Thus, our advection scheme conserves mass up to the rounding errors. In turn, B SO2,a is also constructed with the mass conservation. However, our smoothing procedure leads to slight violation of the mass conservation. 170 We chose to recover this conservation by adjusting B SO2,smo with the scalar adjustment coefficient ν adj : where " global " stands for the area-weighted summation over all model grid cells. .
The ChAP data flow is summarised in Fig. 2.

Parametrisation of chemical sources and sinks
We assume that k in−cl is proportional to cloud fraction c and cloud water path and, in addition, depends on temperature because of the respective dependence of the involved reaction rate constants. As a result, we chose to use where k in−cl,0 , α in−cl , and β in−cl are constants, T 0 = 288 K. In this equation, :: the : dependence on cloud parameters is constructed by assuming that most of the SO 2 oxidation occurs in the cloud-covered part of the model grid cell and taking into account that at the coarse spatial and time scale the cloud water path depends on c approximately as a power function (Eliseev et al., 2013). The :::::::::: dependence :: of ::: the oxidation rate k in−cl dependence on temperature is uncertain as well because this conver-190 sion is not a single-step reaction and depends on solubilities of sulphur substances in water and on the rate of the SO 2 oxidation by peroxide radical. Therefore, it is difficult to relate α in−cl directly to the activation energies of these reactions. However, such activation energies are able to provide an order-of-magnitude estimate for value of this coefficient. For instance, the activation energy value for reaction HSO 3 + H 2 O 2 as listed in Table 1 of (Barth et al., 2000) for typical lower tropospheric temperatures corresponds to α in−cl = 0.05 K −1 . We use this value as a guide below.
Recall that k gas = 0 because of R gas = 0 (Sect. 2.1). Thus, the production of sulphates, R SO4,prod ≡ R in−cl . We will discuss this limitation below (Sect. 6).
We set k SO2,dry and k SO4,dry to constant values (Table 1). Because sulphates wet deposition should depend on precipitation rate p and this dependence is expected to saturate somewhat in a limiting case of very strong precipitation, after some trialand-error procedure we chose where k SO4,wet,0 and p 0 are constants.

Simulations setup
We ran our model for 1850-2000 with the SO 2 emissions data from the CMIP5 (Coupled Models Intercomparison Project, phase 5) 'historical' database (Lamarque et al., 2010) (see also Fig. S3). This database lacks the global gridded data on B SO2 205 but provides the data on B SO4 (Lamarque et al., 2013b). The data on sulphate burden were used to evaluate the performance of our scheme. Both emission and burden data are available as time slices with a step of 10 yr.
The CMIP5 data are recently superseded by the CMIP6 (Coupled Models Intercomparison Project, phase 6) data sets (Hoesly et al., 2018;Turnock et al., 2020). However, because the CMIP5 data are sufficient to validate our scheme, and because we expect that this scheme would need some (but not major) retuning when it is implemented into an Earth System Model, we 210 limit our calculations in the present paper with :: to : the CMIP5 data. We postpone the task to run our scheme with the CMIP6 emissions for the next stage -when our scheme is implemented into EMIC.
In our calculations, we neglect dimethyl sulphide emissions from the ocean, which is an important source of the sulphur dioxide in the marine atmosphere (Warneck, 2000;Surkova, 2002). We do it mostly because they are not available in the CMIP5 forcing data (see https://tntcat.iiasa.ac.at/RcpDb/). Moreover, we neglect other, more minor sulphur sources such as 215 volcanos and the terrestrial biosphere. Thus, we assume that natural sources did not change since year 1850, which is in the CMIP5 protocol is considered as a pre-industrial year. In addition, we note that an implementation of the natural sources for the considered here sulphur compounds, while certainly important, would complicate our scheme. Again, we postponed this task for future work. Therefore, we compared the given year our simulation with the difference of the CMIP5 data for this year from the respective data for year 1850. 220 In addition, we neglected the direct anthropogenic SO 4 emissions into the atmosphere which contribution to the sulphur budget is generally small (Houghton et al., 2001, their Table 5.5). However, a possibility to account for these emissions is already coded in ChAP and may be used in future simulations.
We use the monthly mean ERA-Interim data (Dee et al., 2011) averaged over 1979-2015 to force our scheme. This set up neglects dependence of meteorological variables on time, and, therefore, respective dependencies of spieces advection.

225
In addition, this approach ignores interannual changes of temperature in Eq. (14). However, a similar neglect is embedded into the construction of the CMIP5 SO 4 burdens (Lamarque et al., 2013b). Thus, our approach even makes the evaluation of our scheme more straightforward.

235
To tune our scheme, we follow the procedure which is similar to that used by Eliseev et al. (2013). At first, we tune it manually to achieve a first-guess, reasonable performance. At the next stage, we sample the first 7 parameters listed in Table 1 in the predetermined intervals. The sampling was done by using the Latin hypercube sampling (McKay et al., 1979;Stein, 1987) to insure that the ensemble statistics is unbiased. The sample length is K = 5, 000.

250
Upon completing model runs with each parameter set from this sample, we selected only those runs which fulfil the requirements S k ≥ 0.06, 0.8 ≤ R SO4,prod /D SO2,dry ≤ 1.2 and T SO4 < 7 days. The first requirement is based on the observation that the maximum value of S k is 0.7113, so we choose the simulations which are close to the optimal. The second requirement arises from simulations reported in (Warneck, 2000;Surkova, 2002;Houghton et al., 2001, their Table 5.5), in which sulphur dioxide emissions were almost equipartitioned between production of sulphates and SO 2 deposition. The third requirement is 255 based on typical lifetimes of sulphates in the atmosphere. While it looks redundant taking into account our calculation of s k,m , it is necessary to avoid an overfitting of the observed fields in the regions of small sulphate burdens, which are not so important for climate and ecological applications. Such overfitting in our simualtions tends to bias the model with underestimated SO 4 production (this, in addition, motivated us to implement the requirement on R SO4,prod /D SO2,dry ) and deposition of sulphates, despite of the reasonable SO 4 burden.

260
As a result, 40 simulations were considered as being close to the optimal. The means of parameters over these simulations were considered as a tuned parameter set, and their standard deviations were considered as a measure of uncertainty for these parameters.

Perfomance
The tuned parameter values and their uncertainties are listed in Table 1. Below, only the simulations with the tuned set of 265 parameters is discussed.
We assessed the performance of our tuned scheme by comparing it to the original CMIP5 data ( :::: these :::: data :::: were :::: used :: to :::: tune ::: the ::::::: scheme, ::::: thus, to avoid a circular reasoning, we highlight that it is an evaluation of our tuning procedure rather then of the implemented physics); -ACCMIP phase II simulations (Lamarque et al., 2013a;, which were performed both for the preindus-270 trial and for the present day emissions of aerosols and their precursors to the atmosphere (thus, the difference between these simulations is an analogue to our anthropogenic-only simulations). The caveat, however, is due to difference in prescribed SO 2 emissions between our simulations and the ACCMIP protocol (Table 2).
In addition, we use two datasets based on the assimilation of the available measurements into the chemical-transport models: the  Table 5.5 of IPCC TAR (Houghton et al., 2001). We note that, in contrast to the CMIP5 and ACCMIP datasets, CAMS and EMEP MSC-W were prepared by using the meteorology which changes from year to year making makes such comparison less straightforward. The individual model simulations in the IPCC TAR Table 5.5 were performed in a stationary fashion, with the year-to-year changes in meteorology only due to internal model variability. There-285 fore, we can consider them as also being ran with an 'almost constant' meteorological fields, which makes our comparison with these simulations more straightforward.
Below we first discuss the performance of our scheme with respect to simulation of SO 2 , because it is independent from the SO 4 performance. Then we proceed with a similar discussion of SO 4 simulation which, in contrast, depends on the calculated SO 2 burden. We note that the maximum anthropogenic SO 2 emissions into the atmosphere in the CMIP5 data correspond to 290 the 1980 time slice. However, because the model results for year 1990 are quite similar to those for year 1980, and because the most data are available since year 2000, we use the time slice for year 1990 as a primary model output to compare to the existing data. In such cases, we use the year 2000 time slice as a primary source for comparison.

Simulation of SO 2 burden and near-surface concentration
At the global scale, about half of the SO 2 emissions in our model are consumed by the chemical SO 4 production in the 295 atmosphere, and another half is deposited to the surface in the form of sulphur dioxide (Table 2). This fractions are within the ranges reported in IPCC TAR  Fig. 3a). These burdens are in the lower range of the IPCC TAR Table 5.5-derived values (  Nonetheless, one notes some overestimation of both variables in Europe and some underestimate in south-east Asia. In Europe, the ChAP-simulated q SO2,s also agrees with the EMEP MSC-W data for 2000-2005 (Fig. S6). The larger discrepancy of our simulations in Europe from the CAMS data than from the EMEP MSC-W data is at least partly explained by difference in covered period between the CAMS and EMEP MSC-W datasets. Namely, provided that aerosol emissions in Europe continuously burden for year 1990 is similar to that simulated with the NCAR CCM (Barth et al., 2000), GISS (Chin et al., 2000) and CCCMA (Lohmann et al., 1999) al., 1995) and GISS (Chin et al., 1996) models.

Simulation of SO 4 burden and near-surface concentration
Similar to it was for SO 2 , the total column burden of anthropogenic sulphates monotonically increases until 1980, reaches The magnitude of B SO4 in the North American source regions is basically correct, but during summer the maximum in this region is shifted to the west. The latter feature is not exhibited in winter. The B SO4 magnitudes in the source regions in the Southern Hemisphere are overestimated for the whole year.
Compared to the CAMS reanalysis ( Fig. S4) in major source regions, our model overestimates sulphates burden per unit 360 area in Europe with the larger discrepancy in winter then in summer. The B SO4 pattern in south-east Asian source region is underestimated -this differs to that obtained in the comparison between our ChAP simulation and the CMIP5 database.
The latter difference, at least partly, is due to difference in covering periods (recall, that CAMS is for 2003-2010). Again, the magnitude of the SO 4 burden in the North American source region is realistic in ChAP, but now we see that even the location of maximum is correct. Thus, our previous conclusion about this location is points to some possible shortcomings in the CMIP5 365 dataset. In the Southern Hemisphere, our model overestimates sulphate burdens per unit area in both south African and South American source regions. In addition, B SO4 in the Northern Hemisphere source regions is similar to those reported in the simulations with the NCAR CCM (Barth et al., 2000;Rasch et al., 2000), ECHAM (Feichter et al., 1996;Roelofs et al., 1998), GISS (Chin et al., 2000), and CCCMA (Lohmann et al., 1999) models.
Geographic distribution of near-surface SO 4 concentration, q SO4,s , follows the corresponding distribution of B SO4 (Fig. 8).

370
Identically to that it was for sulphur dioxide, this is a direct consequence of Eq.
(2). In the anthropogenically polluted regions, q SO4,s in the last decades of the 20th century is above 2 µgS m −3 , and in Europe during summer 1990 it is larger than

Simulation of annual SO x deposition 380
Owing to the mass conservation, the global SO x deposition in the model is equal to the applied sulphur dioxide emissions.
Depending on time slice, dry SO x deposition D SOx,dry = D SO2,dry + D SO4,dry explains from 55 to 59% of the total SO x deposition (mostly in the form of SO 2 ), and wet deposition D SOx,wet = D SO4,wet explains another 41 − 45% (only in the SO 4 form by construction ::::: design) ( Table 2, Fig. 9a). The contribution of wet SO x deposition in 1980 and 2000 is also similar to that obtained from ACCMIP (46% and 51%, respectively, Lamarque et al., 2013a) and is within the ranges reported in Table 5.5 385 of(Houghton et al., 2001) for year 1990 (from 37 to 64% with mean of 47% and median of 45%).
Geographic distribution of the total SO x deposition D SOx = D SOx,wet + D SOx,dry is very close to the sulphur dioxide emissions in a given year (cf. Fig. 9a with  . 9b, c and S7a, b). Again, this is a validation for our numerics and for our code rather than for the implemented physics, because total SO x deposition near any source region is controlled by prescribed emissions and by prescribed winds.
More stringent test of the implemented physics is a subdivision of total SO x deposition into wet and dry ones (Figs. 9d-g and S7c-f). This shows that ChAP generally overestimates wet deposition and underestimates dry one relative to the ACCMIP 395 simulations. This is not visible in the global numbers (Table 2, Table 2 in Koch et al., 1999), such increase would change sulphur dioxide and suphate burdens in the troposphere mostly over oceans. We plan to implement this source into our scheme in future.
Another transport-solver related issue is due to implemented smoothing procedure. This implementation is reasoned by neglect of synoptic scales in our transport routine -we use only monthly mean winds. While ChAP is formally linear with 440 respect to horizontal winds (Eq. (5), and, therefore allows averaging over synoptic-scale motions, possible time correlations between synoptic-scale variations of winds and of pollutants burdens make the underlying processes nonlinear (Saltzman, 1978;Branscome, 1983;Petoukhov et al., 2008;Coumou et al., 2011). Further, one may argue that the corresponding mixing length (thus, n smo ) could be made dependent on synoptic-scale kinetic energy (Branscome, 1983;Coumou et al., 2011). At the time being, n smo is a parameter of the scheme, but in future it could become depending on large-scale atmospheric state (e.g.,

445
on the state with time and space scales larger than synoptic ones).
In our scheme we neglected wet deposition of sulphur dioxide. This was done based on the synthesis of simulations listed in Table 5.5 of IPCC TAR (Houghton et al., 2001), in which wet deposition of SO 2 explained no more than 15% of sulphur dioxide budget. Only in CCM1- GRANTOUR (Chuang et al., 1997) and :: in the earlier version of GOCART (Chin et al., 1996) this contribution is from 15 to 20%. While the upper-end values from these papers are not negligible, D SO2,wet is still neglected.

450
Its implementation would probably improve regional performance of our scheme. We acknowledge the neglect of D SO2,wet as a limitation of our scheme.
In addition, is is necessary to highlight that dry deposition of sulphates is still included in our scheme, despite its contribution to the SO 4 budget is also not so important relative to the SO 2 wet deposition at the global scale. For instance, for the IPCC TAR ensemble this contribution is ≤ 25%. The reason for keeping D SO4,wet is mostly numeric: such background deposition 455 avoids a division by zero in regions with very small precipitation rate.
More subtle issue is due to implementation of total cloud fraction in Eq. (14) for in-cloud oxidation rate. Because most sulphur conversion occurs in the lower half of the troposphere, one may argue that low cloud fraction would be a better predictor for this oxidation rate. We tried this option during development of ChAP and found no marked differences (apart somewhat different optimal values of parameters listed in Table 1). Thus, to be in line with the contemporary generation of 485 EMICs, which mostly do not provide cloud fractions for different layers, we kept total cloud fraction in the input to our scheme instead of low cloud fraction.
Another way of improving the calculation of near-surface concentrations of sulphur species is to account for regional differences of H SO2 and H SO4 . This may be done be relating its values to large-scale values of the planetary boundary layer depth (thus, to vertical mixing inside this layer). We postpone this as a possible future extension of our scheme.

490
Performance of the scheme was mainly tested against the CMIP5 data and the ACCMIP simulations. The basic reason for limiting validation to these datasets is due to our neglect of the natural sources of sulphur emissions into the atmosphere.
This limitation is somewhat relaxed in the present paper by comparing to the CAMS and EMEP output and to individual simulations in regions of strong anthropogenic pollution into the atmosphere. Our neglect of natural sulphur sources is also one of the reasons to exclude direct measurements of sulphur burdens (e.g., Aas et al., 2019). Another reason for this exclusion is due to possible complications owing to local features which are likely present in these direct, point-scale measurements.
A meaningful use of such data would need their stratification into background and polluted stations and factoring out such local-scale features, which is beyond the scope of the present study. Exclusion of such point-scale measurements is somewhat ameliorated by their assimilation into CAMS and EMEP. In addition, in future exercises, we plan to replace the CMIP5 forcing by the CMIP6 one (Hoesly et al., 2018;Turnock et al., 2020).

500
A related issue is due to our use of the ERA cloud and precipitation fields rather than those based on direct measurements.
For instance, arguably more reliable cloud fractions may be prescribed from the A-Train satellite observations (Minnis et al., 2011;Frey et al., 2008). Correspondingly, precipitation rate could be derived from the GPCP (Global Precipitation Climatology Project) data (an update from Huffman et al., 2009). Marked differences of the ERA-Interim fields from these data are documented, for instance, by Dolinar et al. (2016), Stengel et al. (2018), and by Nogueira (2020). For instance, the underestimated 505 rainfall rate in Europe in ERA-Interim (Nogueira, 2020) region may be the reason of our overestimate of B SO4 in this region, while the corresponding excessive precipitation rate in the Asian monsoon region could contribute to the underpredicted sulphate burden over south-east Asia. However, we prefer to keep the ERA-Interim cloud fractions and precipitation as forcing fields in our tuning exercise because they are at least dynamically consistent with other forcing fields.
Finally, choice of skill scores is always subjective in the tuning exercises like ours. We reported only one version of skill 510 scores (Eq. (16)). However, we tried another skill scores as well, e.g., either based on the root-mean-square errors (RMSE) in simulations or by limiting the skill scores calculations only to the regions with sufficiently large SO 2 and SO 4 burdens.
The first option (RMSE skill score) provided rather unrobust results. The second options did not resulted to much improved simulation with respect to that reported in Sect. 5. Therefore, we decided to use the skill score as figured in Eq. (16).  sources. In particular, in 1980 and1990, when the global anthropogenic emission of sulphur is at the maximum, global atmospheric burdens of SO 2 and SO 4 account, correspondingly, 0.2 TgS and 0.4 TgS. In our scheme, about half of the emitted sulphur dioxide is deposited to the surface and the rest in oxidised into sulphates. In turn, sulphates are 530 mostly removed from the atmosphere by wet deposition. The lifetime of the sulphur dioxide and sulphates in the atmosphere is close to, respectively, 1 day and 5 days. The differences between our simulations, on one hand, and the CAMS and EMEP MSC-W datasets, on the other, are partly (but likely far from completely) explained by the differences in time intervals covered by our simualtions :::::::::: simulations and by these datasets.

545
-Simplifications in the transport solver.
We plan to relax these limitations during future development of our scheme.
Despite its simplicity, our scheme is able to reproduce gross characteristics of the tropospheric sulphur cycle during the historical period. Thus, it may be successfully used to simulate anthropogenic sulphur pollution in the atmosphere at coarse spatial and time scales. At next stage, we are going to implement it into EMIC and reproduce direct radiative effect of sulphates 550 on climate, their respective indirect (cloud-and precipitation-related) effects, as well as an impact of sulphur compounds on the terrestrial carbon cycle.
Code and data availability. The Fortran code for ChAP as well as all the data used in this paper are available at the ZENODO repository via https://doi.org/10.5281/zenodo.4513909.
by comparing the simulations with the EMEP and CAMS data. A.V.T. prepared the data to run the scheme. All authors contributed to the manuscript revisions after the first draft.
Competing interests. The authors declare that they have no conflict of interest.