MIPAS observations of volcanic sulphate aerosol and sulphur dioxide in the stratosphere

Volcanic eruptions can increase the stratospheric sulphur content by orders of magnitude above the background level and are the most important source of variability of stratospheric sulphur loading. We present a set of vertical profiles of sulphate aerosol volume densities and derived liquid-phase H 2 SO 4 mole-fractions for 2005–2012, retrieved from infrared limb emission measurements by the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) on board of the Environmental Satellite. The MIPAS aerosol dataset has been corrected for a possible altitude-dependent bias by comparison with balloon-borne in situ aerosol measurements at Laramie, Wyoming. The MIPAS data of stratospheric sulphate aerosol is linked to MIPAS observations of sulphur dioxide (SO 2 ) with the help of Chemical Transport Model simulations. We investigate the production of sulphate aerosol and its fate from volcanically emitted SO 2 for two volcanic case studies: the eruptions of Kasatochi in 2008 and Sarychev in 2009, which both occurred in the Northern Hemisphere mid-latitudes during boreal summer. We show that the MIPAS sulphate aerosol and SO 2 data are qualitatively and quantitatively consistent to each other. Further, we demonstrate that the lifetime of SO 2 is well described by its oxidation by hydroxyl radicals. While sedimentation of the sulphate aerosol plays a role, we find that the dominant mechanism controlling the stratospheric lifetime of sulphur after these volcanic eruptions at mid-latitudes is transport in the Brewer-Dobson circulation. Sulphur emitted by the two mid-latitude volcanoes resides mostly north of 30°

composition typically varies between around 70 and 80 %, as obtained by equilibrium calculations (Carslaw et al., 1995) and observations (e.g. Doeringer et al., 2012). E.g. Kleinschmitt et al. (2017), calculating aerosol optical properties, Kremser et al. (2016), calculating sulphur fluxes, and Gao et al. (2007), calculating atmospheric volcanic aerosol loadings, also use a 75 wt% H2SO4-H2O composition in their studies. The imaginary parts of various refractive index datasets in the mid-IR are displayed in Fig. 1. Here the used optical constants by Niedziela et al. (1999) for 75 wt% and 230 K (bold red line) are compared to data at other concentrations and temperatures by Niedziela et al. (1999) (upper panel), and Myhre et al. (2003) (lower panel). This particular data set has been chosen because in an evaluation of optical constants for sulphuric acid, Wagner et al. (2003) have found those datasets to be best consistent with observations in the aerosol chamber AIDA (Aerosol Interactions and Dynamics in the Atmosphere). The spectral range selected for the retrieval (1,216.5-1,219.5 cm -1 ) is situated at the long wavelength end of MIPAS band B as indicated by the two vertical lines in Fig. 1. It lies within one of the atmospheric windows as can be seen by comparison to the limb-transmission curve (light grey) in Fig. 1. We have not chosen the windows at around 830 cm -1 and 950 cm -1 since at 1,220 cm -1 the absorption by H2SO4 droplets is higher and the relative difference between the various sets of refractive indices is lower.  The bold red line indicates the dataset, and the two vertical lines the spectral window used in this study. A simulated limb transmission spectrum for 10 km tangent altitude for standard mid-latitude conditions is additionally plotted in the top row as well as the wavenumbers covered by MIPAS bands A, AB and B.
The retrieval has been set up as a multiparameter nonlinear least-squares fit of the calculated to the measured limb radiances of entire limb-scans (e.g. von Clarmann et al., 2003). Besides the target parameter, namely sulphate aerosol volume densities, further atmospheric fit-parameters of the retrieval are vertical profiles of spectrally interfering trace gases CH4, H2O, O3, and HNO3. While zero initial guess profiles have been used for the volume densities, results from the IMK routine processing are taken for the trace gases . As the atmospheric parameters are represented at denser altitude levels (1 km) than the vertical field-of-view (~3 km) and the vertical tangent point spacing (1.5 km) of MIPAS, constraints on the smoothness of the profile shape are introduced by regularization (Tikhonov, 1963;Steck, 2002).
The retrieval of aerosol volume density is restricted to altitudes up to 33 km and the regularization strength has been adjusted such that its resulting vertical resolution is around 3 to 4 km. To cover instrumental uncertainties a spectral shift parameter and a radiance offset, constant over all wavenumbers and tangent altitudes, has been retrieved simultaneously to the 7 5 10 15 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2017-538 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 16 June 2017 c Author(s) 2017. CC BY 4.0 License. atmospheric quantities. For the analysis in this paper only data at altitudes with averaging kernel diagonal elements larger than 0.05 which are at least 1 km above the lowest tangent height are used.
An overview of the leading error components is presented in Fig. 2, with the assumed parameter uncertainties listed in the caption. The error contributions are estimated from a subset of a few hundred single cases by sensitivity studies using modified parameters or, in case of spectral noise, directly from the retrieval diagnostics. The total error changes with altitude from around 20 % (0.09 µm 3 cm -3 ) at 10 km up to over 40 % (0.005 µm 3 cm -3 ) at 30 km. It is dominated by the uncertainty of the optical constants resulting in 10-20 % error, followed by tangent pointing knowledge with 5-15 %. The error component resulting from spectral noise is rather constant with altitude in absolute terms of volume density and amounts to about 0.01 µm 3 cm -3 . Other instrumental errors that have been investigated but are not listed in Fig. 2 are uncertainties due to the knowledge of the instrumental line shape, and radiometric gain and offset calibration error. In the estimation of the radiative error no offset variations with tangent altitudes were considered, and, thus, compensated for by the retrieval approach.
However, a tangent altitude dependent radiometric offset error caused e.g. by straylight in the instrument cannot be excluded (López-Puertas et al., 2009). We have not handled this uncertainty in the framework of error estimation but we have tried instead to compensate for it through a de-biasing of the dataset based on validation with in situ observations as described in Sect. 3.2.

Figure 2.
Altitude dependent estimated errors for the retrieval of H2SO4-H2O aerosol volume densities. Solid lines indicate the uncertainties used to calculate the 'total' error. Indicated errors are: 'noi': single scan spectral noise, 'T': temperature uncertainty 2 K, 'vmr': 10 % uncertainty of volume mixing ratios of interfering gases, 'htang': 300 m tangent altitude uncertainty, 'Niedziela_240K_60wts': use of optical constants by Niedziela et al. (1999) for 240 K and 69 wt% instead of Niedziela et al. (1999) for 230 K and 75 wt% H2SO4. The dotted curves ('Myhre_243K_65wts') show the results when the optical constants of Myhre et al. (2003) for a temperature of 243 K and a concentration of 65 wt% are used instead of those by Niedziela et al. (1999)  Prior to the retrieval, a deselection of spectra affected by clouds has been performed via application of an established cloud filter method for MIPAS by Spang et al. (2004). To sort out optically thick clouds, but not all aerosol-affected spectra, this cloud filter has been applied with a cloud index limit of 1.7. Due to this loose setting of the cloud-filter artefacts caused e.g. by thin cirrus clouds, polar stratospheric clouds (PSCs) or volcanic ash remain in the dataset, which are all attributed to the retrieved 75 wt% H2SO4-H2O aerosol volume density. Thus, further filtering of affected profiles has been necessary after completing the retrieval.
Two distinct features of strong enhancements with an annual cycle show up in the unfiltered dataset. The first feature is due to strong enhancements in presence of PSCs at the winter pole. To deselect PSC-affected profiles a filter is applied when temperatures in the altitude range 17-23 km drop below a threshold of 195 K polewards of 40°, for the Northern Hemisphere from 15 Nov-15 Apr, and for the Southern Hemisphere from 1 Apr-30 Nov. This temperature represents the nitric acid trihydrate (NAT) existence temperature at around 20 km, under typical stratospheric conditions for nitric acid (HNO3) and water vapour (H2O).
The second feature is assumed to be induced by thin cirrus clouds. It is present mainly in the tropics, at around 25° S-25° N, at altitudes between about 13-21 km. It reaches highest altitudes and is most intense above and in the vicinity of continents, and above the western Pacific. The vertical extent is smallest in boreal summer, and its vertical gradient towards lower aerosol volume densities is relatively strong, with no upward transport being observed. Both, the location in the tropics, in regions of strong vertical motions and convective clouds, and the relatively sharp decrease at higher altitudes towards increasing temperatures suggest that it is connected to the influence of ice particles. The ice-filter for MIPAS data by Griessbach et al. (2016) is applied on all retrieved MIPAS aerosol profiles to reduce the effect of spectra influenced by ice in the present dataset. This method consists of two steps to detect whether MIPAS spectra are influenced by aerosols, ice, clouds, ashes, or a clear sky (Griessbach et al., 2014 and. First aerosols and clouds are identified, using a spectral window region that is sensitive to aerosols and clouds. Then ice clouds and aerosols are discriminated, using spectral windows with contrasting behaviour for ice and aerosols. This is then combined to brightness temperature difference correlations. In our dataset, we consider only retrieved values starting 4 km above the altitude of the uppermost spectrum that was flagged to have been influenced by ice. Further, the ash filter for MIPAS spectra by Griessbach et al. (2014), based on an ash detection threshold function, is applied in the same way as the ice filter, to filter out volcanic ash and mineral dust.

Validation and bias correction
To validate the new dataset, we compare the profiles of MIPAS aerosol volume density to in situ balloon measurements (Deshler et al., 2003). In situ measurements were carried out with laser based aerosol spectrometers (LPCs) from Laramie, Wyoming (41° N/105° W), between 6 and 9am, local time. In Fig. 3, profiles of the balloon measurements are shown. In comparison, MIPAS mean aerosol volume density profiles are presented, selected from a restricted area around Laramie, together with their normalised standard deviations (Fig. 3)   Generally, the aerosol volume densities ( Fig. 3 and 4) are highest in the lower stratosphere and then decrease towards zero at higher altitudes. As the balloon data has a higher vertical resolution, and the retrieval process for MIPAS profiles includes smoothing, the in situ data show finer structures. Compared to the balloon data, the original MIPAS aerosol volume densities show a positive bias in most profiles (Fig. 3) as well as in the mean profile (Fig. 4). This is most easily detectable at higher altitudes where profiles are relatively smooth. The offset amplifies towards lower altitudes (Fig. 4b). Aiming on a reduction of this positive offset, a height dependent de-biasing is performed on all single MIPAS profiles. The de-biasing is based on the in situ measurements carried out with laser based particle counters. MIPAS profiles show a consistent variation with height, compared to the LPC measurements. An additive linear de-biasing is applied, rather than a multiplicative correction, as the offset is expected to be caused by a possible altitude-dependent additive stray light error in the radiances (see Sect. 3.1). The offset-correction is represented in Fig. 4b (dashed line). It is the difference between a linear fit to MIPAS values at 18 to 30 km altitude, and the in situ data (Fig. 4b, pink and blue solid line, respectively). At lower altitudes, where profiles show more variability, this linear fit also suits well. The mean de-biased MIPAS profile ( Fig. 4a) matches the in situ data and lies mostly in the range of the standard error of the mean of the in situ data. Further the absolute and relative differences to the balloon data are reduced significantly ( Fig. 4b and c). Percentage differences are mostly below ± 25 %. For the non-de-biased profile, at altitudes above around 20 km percentage differences increase strongly, due to very low aerosol volume densities, while at lower altitudes percentage differences are below about 100 %.

Time series of MIPAS sulphate aerosol and SO2 for 2005 to 2012
To study the distribution of sulphate aerosol, as measured by  1. In polar regions at altitudes above ~16 km, sulphate aerosol mole-fractions decrease strongly in winter to spring. The pattern is more pronounced in the Southern Hemisphere. This decrease is connected to the polar vortex, where relatively sulphate aerosol free air is transported downwards. Thomason and Poole (1993) reported on very low observed aerosol levels relative to non-vortex air.
2. In both hemispheres, but primarily in the Southern Hemisphere, mole-fractions of liquid-phase H2SO4 are enhanced at around 20-22 km in the mid-latitudes (and partly the tropics) during boreal / aural winter and spring, respectively. satellite measurements are compared to CTM simulations, to study the fate of the emitted sulphur in terms of conversion from SO2 to sulphate aerosol, and its transport and removal at altitudes between 10 and 22 km. As our intention is to study explicitly the sulphur per volcanic eruption, background values per model simulation are set to zero for both SO2 and H2SO4, and no other sources than the volcanically emitted SO2 of one volcanic eruption is included.

Sulphur mass in the Northern Hemisphere mid-and high-latitudes
In this section we aim at testing the agreement between measured SO2 and liquid-phase H2SO4 masses, together with modelled data, in terms of the increase and decline of sulphur emitted by the volcanic eruptions of Kasatochi in August 2008 and Sarychev in June 2009, and the influence of sedimentation radius on the residence time of sulphate aerosol. Good accordance between the modelled and measured SO2 masses is essential to test, by comparison with modelled sulphate aerosol, if the measured aerosol is quantitatively and qualitatively consistent with its measured precursor.
In Table 1, SO2 masses for three altitude regions, as used for the CTM simulations in the present study, are given (labelled 'present study'), together with comparisons to volcanic SO2 masses from the literature. The simulations result in good agreement between measured and modelled SO2. The upper injection limit for the volcanic emissions in the CTM is set to 19 km, an altitude limit derived from comparisons with MIPAS SO2. The main part of SO2 per eruption is emitted into the altitude region from 10 to 18 km, and only few percent of the SO2 masses are injected into altitudes above 18 km. In the case of Kasatochi, our model is run with the SO2 masses by Höpfner et al. (2015), reduced by their given uncertainties. Höpfner et al. (2015) derived volcanic SO2 masses for three altitude regions from 10-14 km, 14-18 km, and 18-22 km by exponential extrapolation of the MIPAS SO2 masses back to the eruption day. This method was applied as in the first month after the eruption MIPAS underestimates the SO2 (Höpfner et al., 2015). Error bars of the extrapolated values were found to be relatively large (Höpfner et al., 2015;presented also in Table 1). In the case of Sarychev, the SO2 mass used in the present study lies below the given error bars.
When comparing the SO2 masses from different studies, it has to be pointed out, that the SO 2 masses are generally not derived for the same altitude regions. Höpfner et al. (2015), Brühl et al. (2015), and the present study are not totally independent from each other, as they are entirely or partly based on the same MIPAS SO2 data by Höpfner et al. (2015). The SO2 masses in our study lie below all studies but Brühl et al. (2015) for Kasatochi, and in the range of the other publications for Sarychev. The wide range of SO2 masses in Table 1 shows the difficulties and uncertainties related to the determination of volcanically emitted SO2.   Höpfner et al. (2015) the given total uncertainty is the sum of the uncertainties per altitude range (Table 3 Fig. 6a-c, the sulphur mass is shown for SO2 and sulphate aerosol separately, while Fig. 6d-f is the total sulphur contained in SO2 and sulphate aerosol. The total simulated sulphur mass is not influenced by chemical sulphur removal, but only by removal due to transport by advection and sedimentation. To analyse the measured and simulated data, datasets of sulphur mass densities (SMD = mass per unit volume) are resampled on a common grid with 1 km vertical spacing and a horizontal resolution that equals the model grid. On this new grid, the same data basis is used for the measured and simulated data, dismissing all 'grid cells' for which either only MIPAS or only CTM data are available. For MIPAS aerosols, SMDs are calculated from the primarily retrieved volume densities, using an assumed aerosol density of 1,700 kg m -3 , and a binary solution of 75 wt% H2SO4-H2O, while for MIPAS SO2 and the CTM data SMDs are calculated from the measured and simulated mole-fractions. Sulphur masses are then derived from 5-days running zonal means of SMDs, by multiplication with the corresponding air volume of the new grid.
Generally, when calculating an integrated mass, high data coverage is crucial to prevent underestimation, therefore we use a 5-days running zonal means. Zonal mean values, used to calculate sulphur masses, are derived using a method of increasing area averaging (see Appendix), to reduce the bias of mean values due to a non-uniform data coverage. Even though high data coverage is very important, we use a method that dismisses available data and information, as the same 5 10 15 20 basis of available values is used for MIPAS and the CTM. This is appropriate when analysing the agreement between the data. Data is especially dismissed for the CTM. Thus we also provide some information on modelled sulphur masses derived from the non-co-located data (Fig. 6a-c). The impact of missing data is strongest in the lowermost altitude region presented here. For MIPAS this is mainly due to the presence of clouds and ash, which were filtered out using the cloud filter by Spang et al. (2004) in the case of SO2 and partly filtered out in the case of aerosol, and additionally the ice and ash filters by Griessbach et al. (2016 and2014, respectively) for the aerosol retrieval. The CTM has low data coverage at lower altitudes due to its isentropic vertical grid. Interpolation to geometric heights starting at 10 km, produces missing values at altitudes up to 13 km, especially in mid-to high-latitudes, and in the presence of strong vortices.
To ease visual comparisons of measured and modelled sulphur mass in Fig. 6, a constant background is added to the model results, as only volcanic sulphur is considered in these simulations. The background mass is chosen considering the mass derived by MIPAS before the volcanic eruption in the region of interest, per altitude and latitude bin. This does not necessarily represent normal background conditions, but unmasks the anomalies caused by the volcanoes.
Concerning the measured and modelled SO2 masses after the eruptions of Kasatochi and Sarychev (Fig. 6), comparisons show that until about one month after the eruptions, the SO2 mass is by far underestimated by MIPAS. This underestimation of SO2 was stressed by Höpfner et al. (2015), when comparing MIPAS SO2 to measurements by the Microwave Limb Sounder (MLS), on board Aura . It is mainly due to the presence of particles, that hinders MIPAS SO2 measurements in largely eruption-affected air-parcels and causes a sampling bias towards less volcano-affected air parcels. Through our model simulations we confirm this bias, and the related time scale found by Höpfner et al. (2015).
The decay of SO2 is well simulated by the CTM, in comparison to the MIPAS measurements. Only oxidation by OH is considered in the model, and we see that the decay of SO2 can adequately be described by this mechanism. Other processes, as decay by photolysis or reaction with O are not considered, and following the good agreement between measurements and model results, can be neglected at the temporal scale and region of interest. Inside volcanic plumes chemistry interactions might lead to changes in SO2-lifetimes (Bekki, 1995). When a high amount of SO2 gets depleted by hydroxyl radicals, the concentration of the radicals might decrease, which could reduce the speed of further depletion. The good accordance between MIPAS measurements and CTM simulations, which do not account for any feedback on the OH concentrations, suggests, that even if such interactions occurred, they did not produce a strong impact on the timescale of months and larger spatial scales.
To investigate the effect of particle sedimentation on the residence time of sulphur after the volcanic eruptions, model simulations with different effective sedimentation radii are performed, including one simulation without any sedimentation.
The radii lie in the range of aerosol size distributions, as observed by Deshler et al. (2003) and Deshler (2008) for volcanically perturbed periods, and one constant radius is applied for all H2SO4 droplets per simulation. Fig. 6 shows the influence of varying the gravitational settling between no settling, and radii of 0.1, 0.5, and 1 µm. The amount of sulphate aerosol removed by sedimentation increases with growing particle size, while the time needed for the removal increases for smaller sedimentation radii. The sulphur mass contained in liquid-phase H2SO4 from a simulation with a particle radius of to the decrease of measured aerosol. Here, sulphate aerosol simulated with a radius of 1 µm compares better. A larger sedimentation radius seems more appropriate at lower altitudes, as heavier particles can settle faster, and can be removed more rapidly than smaller and lighter particles. These can float in the atmosphere or undergo ascent. The particle size distributions of aerosols can further show natural variation for different volcanic eruptions; therefore some differences in the agreement between modelled and measured data when studying different volcanic eruptions can be expected. In general, we conclude that a particle radius of 0.5 µm gives a satisfactory fit between the measurements and simulations for the purpose of studying sulphur mass and sulphur transport in the Northern Hemisphere. Hence, we base all following model results on the CTM runs with a sedimentation radius of 0.5 µm.
We conclude from the comparisons between measured and simulated SO2 and sulphate aerosol, that the amplitude of the peak of liquid-phase H2SO4 and its removal from the studied altitude regions, as measured by MIPAS, is consistent with the measured SO2, both qualitatively and quantitatively. In the model, sulphur is released from SO2, due to its reaction with OH, and sulphate aerosol is consequently built. Modelled SO2, that fits well to MIPAS SO2 measurements in terms of amplitude and decay, releases sulphur and builds H2SO4 that in turn matches well to MIPAS sulphate aerosol in terms of amplitude and decrease.
Further, we find that the dominating process on the fate of volcanic sulphur is transport by the Brewer-Dobson circulation out of the region of interest. This becomes obvious when comparing the total modelled sulphur with and without sedimentation to the observed sulphur mass (Fig. 6d-f). In the case of the CTM this excludes all influence by chemical reactions on the removal of volcanic sulphur. Even though consideration of sedimentation of sulphate aerosol with a sedimentation radius between about 0.5 and 1 µm further improves the agreement between model results and observations in 10-22 km altitude, the modelled sulphur mass without sedimentation already compares rather well with the measured sulphur mass. Hereby we see that the removal is dominated by advection rather than sedimentation.
In the altitude region of interest, from around 10 to 22 km height, supplementary processes, as the photolysis of gasphase H2SO4, that is important at altitudes above 30 km (Vaida et al., 2003;Brühl et al., 2015), or a meteoritic dust sink (Brühl et al., 2015), are not considered. Other processes, such as the evolution of sulphate aerosol through microphysical processes, as nucleation, coagulation, condensation, or evaporation, and sedimentation of particles with different sizes can play a role in our region of interest. However, comparisons of simulations and measurements show that these processes are not essential to study the fate of sulphur emitted by Kasatochi in 2008 and Sarychev in 2009.

Sulphur transport
The Kasatochi eruption directly injected a large amount of SO2 especially into the altitude region between 10 and 14 km (Table 1)    In Fig. 9 and 10 time series of latitudinally resolved mole-fractions show the transport of SO2 and sulphate aerosol at different altitudes. We present time series of 5d running zonal mean mole-fractions for the Northern Hemisphere, at altitudes from 10 to 22 km, both for MIPAS measurements and CTM simulations. As the model has very poor coverage at 10 km, results are not shown for the model at this altitude. For the eruption of Kasatochi, the separation of the plume and downward transport of the upper part is notable at mid-to high-latitudes, most easily visible for sulphate aerosol (Fig. 10).   Figure 10. As Fig. 9, but for sulphate aerosol. (right) Black to white contour-lines denote 5, 25, 50, 100, 150, and 200 pptv. Both, in the measurements and simulations, most of the sulphur contained in SO2 and sulphate aerosol stays north of 30° N at lower altitudes up to around 16 km ( Fig. 9 and 10). Especially at low altitudes we find a mixing barrier at ~30° N, with a strong gradient between low values in the tropics and high values in the extra-tropics, which weakens towards higher altitudes. This gradient is due to the subtropical jet stream, and is most easily detectable in the contour lines shown for modelled liquid-phase H2SO4 (Fig. 10, right panel), but similar patterns are observed by MIPAS. At around 16-18 km, especially in the longer lived sulphate aerosol, this forms a "tongue" of relatively high mole-fractions, which persists over a longer period than in the surrounding latitudes. An additional transport process starts at an altitude of about 18 km in the case of Kasatochi and ~16 km in the case of Sarychev ( Fig. 9 and 10). At these altitudes, both the MIPAS measurements and CTM simulations show that sulphur from mid-latitude volcanic eruptions can reach the tropics, predominantly in the form of sulphate aerosol. In the tropics sulphur is then lifted in the 'tropical pipe' (denotation following e.g. Plumb 1996), and can reach the stratospheric 'overworld' (denotation following e.g. Hoskins, 1991), also seen in Fig. 8. pattern. They find, that at 360-400 K potential temperature, the southward transport was primarily caused by anticyclonic Rossby wave breaking, intensified by the Asian summer monsoon during Northern Hemisphere summer. Above 400 K less aerosol is transported into the tropics. They further show an 'aerosol hole' in the anticyclone, surrounded by aerosol-rich air.
Following Wu et al. (2017), a strong subptropical jet in combination with weak Rossby wave breaking events would hinder the southward transport of the volcanic plume during winter conditions. The weaker southward transport in the case of the Kasatochi eruption that starts at higher altitudes, compared to the Sarychev eruption, could be due to the eruption having been later during the monsoon season, leading to enhanced southward transport by the Asian summer monsoon for a shorter period of time. In model studies of the Sarychev eruption, Haywood et al. (2010) find that sulphate aerosol is transported around the entire globe in around 14 days. We see that the bulk of the aerosol that moves southwards reaches the equator about 2 to 3 months after the volcanic eruption. To some small extent this sulphur crosses the equator and can hereby influence the sulphur content of the Southern Hemisphere (see also Wu et al., 2017).