Effects of land use and anthropogenic aerosol emissions in the Roman Empire

. As one of the ﬁrst transcontinental polities that led to widespread anthropogenic modiﬁcation of the environment, the inﬂuence of the Roman Empire on European climate has been studied for more than 20 years. Recent advances in our understanding of past land use and aerosol-climate interactions make it valuable to revisit the way humans may have affected the climate of the Roman Era. Here we estimate the effect of humans on some climate variables in the Roman Empire at its apogee, focusing on the impact of anthropogenic land cover as well as aerosol emissions. For this we combined existing land 5 use scenarios with novel estimates (low, medium, high) of aerosol emissions from fuel combustion and burning of agricultural land. Aerosol emissions from agricultural burning were greater than those from fuel consumption, but on the same order of magnitude. Using the global aerosol-enabled climate model ECHAM-HAM-SALSA, we conducted simulations with ﬁxed sea-surface temperatures to gain a ﬁrst impression about the possible climate impact of anthropogenic land cover and aerosols in the Roman 10 Empire. While land use effects induced a regional warming for one of the reconstructions caused by decreases in turbulent ﬂux, aerosol emissions enhanced the cooling effect of clouds and thus led to a cooling in the Roman Empire. Quantifying the anthropogenic inﬂuence on climate is however challenging since our model likely overestimates aerosol-effective radiative forcing and prescribes the sea-surface temperatures.

on climate already in AD 1850. But when did anthropogenic aerosol emissions start to change the climate? Is it possible that locally significant changes occurred already thousands of years ago?
In the first century AD, the Roman Empire was at its apogee, with a "surprisingly high standard of living" for antiquity (Temin, 2006). Around the same time, the Han dynasty in China also controlled a large part of the world's population. These two empires had comparable spatial extents and population sizes (Bielenstein, 1986;Scheidel, 2009). Although the global 5 population was approximately 2.6 and 5.5 times smaller in AD 1 than in AD 1700 and AD 1850, respectively (Klein Goldewijk et al., 2017), industrial activites between 100 BC and AD 200 in both Europe and East Asia already left imprints in ice cores and sedimentary records, for example heavy metals and 13 C-enriched methane (Hong et al., 1996;Brännvall et al., 1999;Sapart et al., 2012). In cities and towns during antiquity, smoke emissions were sufficiently severe to darken buildings and to enforce laws against air pollution (Makra, 2015). 10 Of course, these activities are by no means comparable to the anthropogenic impact on climate under present-day conditions. Nevertheless, we hypothesise that these anthropogenic activities could already have had an influence on climate on a continental scale. The goal of this study is to gain a first impression on whether or not anthropogenic land cover change and aerosol emissions could have influenced climate already in antiquity. We use the global aerosol-climate model ECHAM-HAM-SALSA for this assessment. Since the sea-surface temperatures (SST) are prescribed in all simulations, ocean-atmosphere feedbacks are 15 disabled and the temperature and precipitation responses are dampened. Nevertheless, the results provide valuable information about changes in variables such as surface albedo, turbulent flux, ARE or CRE.
To analyse the influence of land use on climate, we compare a control simulation without land use to simulations including crop and pasture areas representative for AD 100 (Sect. 3.1). Two different reconstructions of anthropogenic land cover were used. Since ECHAM-HAM-SALSA does not calculate a full carbon cycle, we only investigate biogeophysical effects and 20 not the biogeochemical effects. The impact of anthropogenic land cover on secondary organic aerosol (SOA) precursors is simulated.
To assess the impact of anthropogenic aerosol emissions, we first constructed three possible scenarios (called "low", "intermediate", and "high", respectively) of emissions representative for AD 100. We considered emissions from fuel consumption, crop residue burning, and pasture burning. The magnitude of these new aerosol emissions is compared with that of natural 25 fire emissions in Sect. 3.2 and with that of anthropogenic emissions in AD 1850 in Sect. 3.3. Simulations with and without anthropogenic aerosol emissions were compared to quantify their impact on certain climate variables (Sect. 3.4). the political boundaries of the Roman Empire at this time also occurred outside of it. Since the main focus of our study lies on the anthropogenic influence on climate, and humans near but outside the Roman Empire also carried out agriculture and emitted aerosols, our main conclusions should not be affected by the exact geographical definition. Similarly, some boundary conditions (vegetation, fire emissions) refer to a somewhat earlier period than AD 100 (e.g. AD 1) due to data availability, but concerning the large temporal uncertainties when going so far back in time, we do not consider this to be an issue. 5 Associated with these experiments, we needed data for i) the boundary conditions (Sect. 2.2), ii) anthropogenic land cover change (Sect. 2.3), iii) natural aerosol emissions (Sect. 2.4), and iv) anthropogenic aerosol emissions (Sect. 2.5, 2.6, 2.7).

Model
To study potential anthropogenic effects of land cover and aerosols, we used the global aerosol-climate model ECHAM6.3-HAM2.3-SALSA2.0. SALSA stands for Sectional Aerosol module for Large Scale Applications (Kokkola et al., 2008;Bergman 10 et al., 2012;. The aerosol size distribution is described with 10 size sections. The aerosol species black carbon (BC), organic matter (OM = 1.4 · OC, where OC stands for organic carbon), sulfate (SO 4 ), dust, and sea salt are considered. Particles below r = 25 nm comprise exclusively of OM and/or SO 4 , whereas larger particles can contain any aerosol species.
The land component of ECHAM-HAM-SALSA is called JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in 15 Hamburg; Raddatz et al., 2007). In our simulations, heterogeneity in each grid box is represented by geographically varying fractions of 12 different plant functional types. The implementation of anthropogenic land cover change in JSBACH is described in Reick et al. (2013).

Boundary conditions
The greenhouse gas concentrations follow Meinshausen et al. (2017). Both the greenhouse gas concentrations and the orbital 20 parameters were averaged over AD 50-150 (Table S1). To prescribe natural vegetation fractions (Sect. 2.3), SST, and sea ice concentrations (SIC) we used output from a simulation with the Earth System Model MPI-ESM (Bader et al., 2019, in review, Fig. 1), called MPI_no_LCC in the following. The annual cycle of SST and SIC was derived by averaging the MPI_no_LCC output from AD 50 to 150. The MPI-ESM model has the same atmospheric core (ECHAM) and uses the same vegetation model (JSBACH) as ECHAM-HAM-SALSA. MPI_no_LCC ran from 6000 BC to AD 1850 and considered slow forcings, i.e. 25 changes in greenhouse gases (CO 2 , CH 4 , N 2 O) and orbital parameters, but no anthropogenic land cover change. Vegetation was calculated dynamically (Brovkin et al., 2009).
For the simulated sulfur cycle and SOA calculation, climatologies for the following oxidants are needed: H 2 O 2 (only for sulfur), O 3 , OH, and NO 3 . In general, it is uncertain how these oxidants have changed over the past millennia. OH concentrations seem to be relatively stable (Pinto and Khalil, 1991;Lelieveld et al., 2002;Murray et al., 2014). In contrast, ice core measure-30 ments suggest that H 2 O 2 has increased by > 50% over the last 200 years (Sigg and Neftel, 1991;Anklin and Bales, 1997).
Modelling studies suggest that also NO 3 and O 3 have increased since the pre-industrial (Khan et al., 2015;Murray et al., 2014). For O 3 , Pinto and Khalil (1991); Crutzen and Brühl (1993) found that the changes are relatively small between the glacial and pre-industrial; however, the ozone profile shows larger changes between pre-industrial and present-day conditions (Crutzen and Brühl, 1993). In line with this, we expect that differences in oxidant concentrations are larger between AD 1850 and present-day than between AD 100 and AD 1850 due to large anthropogenic emissions of various gases since AD 1850. Therefore, we used climatological monthly mean mixing ratios of oxidants representative for AD 1850 conditions (Fig. 1).
They were derived from simulations with the Community Earth System Model version 2.0 (CESM2.0) Whole Atmosphere 5 Community Climate Model (WACCM 1 ).

Vegetation and land cover change
In our simulations with ECHAM-HAM-SALSA, natural vegetation was not dynamic. The coverages of different natural vegetation types representative for AD 100 were taken from MPI_no_LCC (Fig. 1). These natural vegetation fractions are fixed over time. They are from an earlier year than AD 100 (end of year 10 BC) because vegetation around AD 1 was used to calculate For the sensitivity simulations in ECHAM-HAM-SALSA, we used two different reconstructions to estimate anthropogenic land cover fractions around AD 100: the anthropogenic land cover reconstructions from KK11 (an update of KK10; Kaplan et al., 2011Kaplan et al., , 2012, and the reconstructions from the HYDE database version 3.1 (Klein Goldewijk et al., 2010 in the following). The empirical model of KK11 assumes that per capita land use declines over time. In contrast, the reconstruc-15 tion of cropland and pasture from HYDE11 assumes a nearly constant per capita land use hindcasting approach with allocation algorithms that change over time (Klein Goldewijk et al., 2017). As a consequence, the anthropogenic land cover fraction in the past is considerably higher in the estimate by KK11 compared to the estimate of HYDE11.
KK11 provides information about the fraction of a gridbox subject to anthropogenic land use. For our simulations, we interpolated the values for AD 1 from a 0.5°× 0.5°grid to the Gaussian grid of ECHAM-HAM-SALSA and assumed that 20 the same fraction of natural vegetation is converted to crop and pasture in equal shares. This seems to be a good first order approximation: in the reconstruction of HYDE11, which estimated pasture and crop areas separately, they contribute each roughly 50%, both when averaged over the whole world and when averaged over the study domain. However, there are of course regions (both in HYDE11 and in reality) where either crop or pasture dominated.
HYDE11 provides information about the area of crop and pasture (km 2 per gridcell) on a 5' grid. We divided these areas 25 by the maximum land area available per gridcell (also from HYDE11) to get fractions of crop and pasture per land area and interpolated the values to the Gaussian grid of ECHAM-HAM-SALSA.
The anthropogenic land cover fractions were scaled to the fraction of the vegetable part of the model grid. Due to inconsistencies between the data sets and the model setup, including differences in the land sea mask, the actually applied land use changes are smaller than prescribed in the original data sets. In our model, the total crop areas in the study domain amount to 30 5.46 · 10 5 km 2 and 13.8 · 10 5 km 2 for HYDE11 and KK11, respectively (i.e. 47% and 21% lower than the original estimate).
For pasture, the total areas add up to 4.75 · 10 5 km 2 and 13.8 · 10 5 km 2 , respectively (47% and 12% lower than the original estimate). For comparison, the total land area in the study domain is 163 · 10 5 km 2 in our model. Part of the underestimation in both datasets is due to the binary land sea mask of JSBACH, which can lead to an underestimation along the coast lines: no crop or pasture can grow in gridboxes that are considered to be ocean in JSBACH, even if the anthropogenic land cover fraction from the original estimate would be larger than zero at this location. The underestimation is considerably more pronounced for KK11 since some areas that are subject to land use in this reconstruction are hardly hospitable to plants in JSBACH. Examples 5 are the Arabian Peninsula, which is and presumably was in reality mainly covered by desert, and parts of North Africa. For instance, the area subject to land use between 35°E and 50°E, 20°N and 35°N (roughly corresponding to the part of the Arabian Peninsula lying within our study domain) amounts to 12.5 · 10 5 km 2 in the original estimate of KK11, but only to 2.28 · 10 5 km 2 in the model. To estimate SOA formation from biogenic sources, a volatility basis set (VBS) was used (Kühn et al., 2019, in preparation;Mielonen et al., 2018;Stadtler et al., 2018). The VBS uses three volatility bins following Farina et al. (2010), which are clas-15 sified as non-volatile (log(C*) = -inf), low-volatile (log(C*) = 0) and semi-volatile (log(C*) = 1), where C* is the equilibrium vapour concentration in g m −3 over a flat surface at 298 K. Gas-to-particle partitioning is computed using a kinetic algorithm after Jacobson (2002). The volatile SOA precursors (monoterpenes and isoprene) are converted into SOA forming species using a simplified one-step oxidation chemistry (Sect. 2.2) and grouped into the volatility bins using pre-defined partitioning coefficients. The emission strengths of the biogenic SOA precursors were calculated online based on the vegetation in the 20 model (Henrot et al., 2017). As a consequence, all ECHAM-HAM-SALSA simulations that include anthropogenic land use account for the effect of these land cover changes on SOA precursor emissions (Table S2). DMS emissions from terrestrial sources were set to present-day values (Pham et al., 1995); the emissions are very low compared to oceanic DMS emissions, though.

Natural aerosol emissions
Fires have played an important role in shaping the composition and structure of Mediterranean vegetation communities 25 (Naveh, 1975). To simulate past natural fire emissions, we used a stand-alone version of the carbon and vegetation dynamics sub-model of JSBACH (CBALONE;called CBALANCE in Wilkenskjeld et al., 2014) together with the fire submodel SPIT-FIRE (Thonicke et al., 2010;Lasslop et al., 2014;Rabin et al., 2017). CBALONE-SPITFIRE needs forcing data related to the vegetation and the atmosphere at a daily time resolution as well as a starting point for the carbon pools. To drive CBALONE-SPITFIRE we used data from MPI_no_LCC (illustrated in Fig. 1) from which output at a high temporal resolution was saved 30 for a few selected 30-year periods (among them a period around AD 1 and one around AD 1835). The driving data for the AD 100 CBALONE-SPITFIRE simulations thus represents a slightly earlier period (around AD 1). However, we do not expect a significant difference between AD 1 and AD 100 since the forcing in MPI_no_LCC is very similar for the two periods. For simplicity we thus refer to the AD 1 fire emissions as those of AD 100. The 30-year period around AD 1 was repeatedly used to drive both the spin-up (≈ 100 years) and the analysed simulated period of CBALONE-SPITFIRE. The emissions have a daily resolution and show interannual variability.
We also calculated fire emissions around AD 1835 (not shown). These emissions were not used as input for our Roman Empire ECHAM-HAM-SALSA simulations but for comparison with the fire emissions from van Marle et al. (2017, Sect. S9). We again used high time resolution output from an MPI-ESM simulation (hereinafter MPI_LCC) to drive CBALONE-SPITFIRE. 5 The only difference between MPI_no_LCC and MPI_LCC is that the latter considers anthropogenic land cover change.
Next to carbon pools and variables related to the vegetation or the atmosphere, the fire model depends on the following additional inputs: i) lightning; ii) population density; iii) a regionally varying anthropogenic influence factor a n .
To estimate the lightning frequency, we tested the parameterisation by Magi (2015) that is based on convective precipitation fluxes. We found that the such-derived lightning frequency is very similar for AD 100 and AD 1835. However, the simulated 10 lightning frequency differs from the standard lighting frequency used by the MPI-ESM implementation of SPITFIRE (e.g. different spatial pattern; not shown), which is based on the Lightning Imaging Sensor/Optical Transient Detector (LIS/OTD) as described in Lasslop et al. (2014). The difference between the simulated and observation-based lightning frequency is probably associated with the high uncertainties in the parameterised convection scheme. Furthermore, the parameterisation by Magi (2015) works better for convective mass flux through the 0.44 hybrid-sigma pressure level than for convective precipitation, 15 but the former was not available from the MPI-ESM simulations. Thus, we decided to use the observationally derived lightning frequencies by LIS/OTD for both AD 100 and AD 1835 ( Fig. 1) in order to not introduce a large bias in natural fire ignitions.
For the CBALONE-SPITFIRE simulations around AD 1835, we averaged the population density from HYDE11 over the years 1830 and 1840. Furthermore, we considered anthropogenic land cover change when calculating the fire emissions; the same changes were used as in MPI_LCC. For AD 100, we used CBALONE-SPITFIRE to calculate natural fire emis-20 sions ( Fig. 1), assuming no anthropogenic influence. Hence, no anthropogenic land cover change was considered in these CBALONE-SPITIRE simulations. Furthermore, the population density was set to 0. The ECHAM-HAM-SALSA simulations LCC_HYDE_low, LCC_HYDE_int, and LCC_KK_high (Sect. 2.8) account for anthropogenic aerosol emissions, which include agricultural burning. Using the same natural fire emissions in these simulations as in the other simulations (no_human, LCC_HYDE, and LCC_KK) would lead to an overestimation in total aerosol emissions since natural aerosol emissions 25 should not occur where now crop or pasture grows. Therefore, we reduced the natural fire emissions of the simulations LCC_HYDE_low, LCC_HYDE_int, and LCC_KK_high offline to account for regions subject to anthropogenic land use (Table S2; Fig. 1). As a first order approximation, the natural fire emissions calculated with CBALONE-SPITFIRE were multiplied with (1 − l), where l is the fraction of the vegetated area per gridbox that is covered by crop and pasture.
Humans impact fires directly by ignitions and fire fighting as well as indirectly, e.g. by forest management and landscape 30 fragmentation (Arora and Melton, 2018). It is very likely that the relationship between humans and fires has changed in the past centuries and millennia as a consequence of large cultural and political changes. This is neglected in SPITFIRE since a n , which reflects cultural differences, does not change over time. Instead, a n (and thus the number of anthropogenic ignitions) in SPITFIRE is tuned towards present-day fires where satellite data is available. For estimating the fire emissions in AD 1835, we nevertheless used CBALONE-SPITFIRE with the standard a n , as it was done by the FIREMIP project for calculating fire 35 emissions from 1750 to today (Rabin et al., 2017). For calculating the (natural) fire emissions around AD 100, a n is irrelevant because the population density is set to zero.

Aerosol emissions from fuel consumption
Many variables influencing the anthropogenic aerosol emissions are highly uncertain. Therefore, three sets of variables respectively leading to low, intermediate and high aerosol emissions were estimated for the Roman Empire based on literature. Figure   5 3 illustrates how much these scenarios differ.
Anthropogenic emissions associated with both fuel consumption and agricultural burning were likely to have had pronounced regional variations. Despite this variability, we tried to estimate values representative for the whole of our study domain for variables such as fuel consumption or fuel load.
We treated the anthropogenic emissions in the same way as natural fire emissions, except for the emission height. In contrast 10 to the natural fire emissions, for which the simulated emission profile depends on the planetary boundary layer (Veira et al., 2015), we emitted the anthropogenic aerosols at the surface.
For fuel consumption, the aerosol emissions EM of species i [kg aerosol m −2 s −1 ] were estimated using the following equation: 15 where cons [kg dry_fuel capita −1 s −1 ] is the fuel consumption per capita, popd [capita m −2 ] is the population density, and EF is the aerosol emission factor of species i [kg aerosol kg −1 dry_fuel ]. In the following, we derive estimates for these three variables for the different emission scenarios. We assume that the three variables are independent when calculating the emissions.

Fuel consumption per capita
In the Roman Empire, people produced aerosol particles by burning several types of fuel for different purposes such as cooking, 20 residential heating, heating bath houses, iron production, glass making, pottery production, or cremation (Malanima, 2013;Veal, 2017;Mietz, 2016;Janssen et al., 2017). For iron production, high temperatures are needed, which were only achieved by burning charcoal (Janssen et al., 2017). For other purposes, also wood or agricultural waste products (e.g. olive pits or dung) were used (Mietz, 2016). The Roman Empire consisted of different regions (e.g. rural versus urban; wetter versus drier climate) with different fuel consumptions per capita and different fuel strategies. Presently, it is therefore not possible to estimate the 25 fuel consumption with large confidence. Malanima (2013) estimates that 1-2 kg of wood were consumed per capita per day in the ancient Mediterranean region.
His numbers refer to fuel for residential heating and cooking, while he states that industrial contributions were negligible. In present-day developing countries that mainly depend on fuel wood for producing energy (sometimes used as surrogates for past conditions), typical estimates of domestic fuelwood consumption also range from 1 to 2 kg per capita per day according to  (2010,2011) and median of fuelwood consumption over the countries considered are approximately 1.5 kg per capita per day (assuming 15% moisture content).
Based on the quantitative model from Pompeii, Veal (2017) applied two extreme scenarios for Rome, where the low and the high estimates for fuel consumption are 1 and 2 t per capita per year. This corresponds to 2.7 and 5.4 kg per capita per day, respectively, which is 2.7 times larger than the estimates by Malanima (2013). The fuel estimates by Veal (2017) however 5 refer to fuel (i.e. the sum of wood and charcoal) in contrast to the estimates by Malanima (2013), which refer to wood (either burnt directly or used for charcoal making; he assumes little contribution from the latter). This distinction is important when calculating emissions because several kilogrammes of wood are needed to make one kilogramme of charcoal, making the difference between the two estimates even larger. While Veal's model derived for Pompeii might give reasonable results for Rome, it is likely not applicable to all parts of the Roman Empire, e.g. in the countryside. Nevertheless, it shows that the fuel 10 consumption in the Roman Empire might have been substantially larger than suggested by Malanima (2013), since the estimates of Veal (2017) account for all types of fuel consumption, including e.g. baths and industrial activites. Recently, Janssen et al.
(2017) calculated the wood consumption for the city Sagalassos (2500-3500 inhabitants) to range between 0.6 and 0.8 kg per capita per day for local pottery production and 1.3-3.4 kg per capita per day for heating the bath (oven dry wood). Although Sagalassos might differ from other places, this indicates that the neglect of non-residential sources by Malanima (2013) might

Population size
We base the population density of our study on HYDE11 for the year AD 100 (for sources see Klein Goldewijk et al., 2010. We divided the population counts from HYDE11 (5' grid) by the gridbox area before interpolating the such derived population densities to the Gaussian grid of ECHAM-HAM-SALSA. Using this approach, the HYDE11 population size between 10°E and 50°E, 20°N and 60°N is around 55 million people (the number of people living within the political boundaries of 5 the Roman Empire would be somewhat smaller). However, the population size of the Roman Empire is still debatable since there is disagreement what the census tallies represent. The number of HYDE11 lies in the range of the so-called "low count" scenario (Scheidel, 2008), but there are also proponents of a "middle count" and a "high count" hypothesis (Hin, 2013;Scheidel, 2008). With the "middle count" approach, Hin (2013) arrived at 6.7 million (range between 5 and 10 millions) free citizens in Italy compared to approximately 4 million with the "low count" for 28 BCE. Assuming a constant ratio between free citizens 10 and all people and a similar scaling factor outside Italy, the population in the Roman Empire would be a factor of 1. 25-2.5 higher in the "middle count" approach compared to the "low count" approach. The "high count" would result in a roughly 3 times larger population size, i.e. more than 100 million people (maybe up to 160 million) living in the whole Empire if we assume that the population densities in other parts of the Empire were similar to Italy (Scheidel, 2008(Scheidel, , 2009.
Based on these different literature values, we used the estimate from HYDE11 for our low emission scenario. For the 15 intermediate and the high emission scenarios, we decided to multiply the population densities of the HYDE database with factors of 1.5 and 2.5, respectively.

Aerosol emission factors
Last but not least, we needed to estimate aerosol emission factors from burning biofuel. We compiled an overview of BC, OC, and SO 2 emission measurements from different studies (Sect. S10, Table S17). Composite estimates were not considered. We 20 treated BC and elemental carbon (EC) to be the same. We grouped the measurements according to the fuel type, i.e. wood (key 1 in Table S17), agricultural waste (key 2), charcoal burning (key 5), and charcoal production (key 6). Woody agricultural waste was counted as wood. We neglected coal as a fuel type although it was widespread in Roman Britain (Smith, 1997). This is justified since the centres of the classical civilisations (especially the Mediterranean region) were not rich in coal (Malanima, 2013). A few measurements refer to PM10 (i.e. aerosol/particulate matter with aerodynamic diameters < 10 µm) 25 or PM4 instead of PM2.5, but the difference is usually small (a few percent between PM10 and PM2.5 in Turn et al., 1997).
Similarly, we neglected the difference between SO 2 and SO x since the latter is dominated by SO 2 .
For the different sectors (i.e. wood, agricultural waste, charcoal burning, and charcoal production), we calculated a lower, an intermediate, and an upper estimate for EF after the following procedure: -For all measurements, the mean, the standard deviation, and the number of samples N were collected.

30
-If N was larger than 1, but the standard deviation was not given and could not be calculated (e.g. from plots, confidence intervals, or data in the Supplementary Material), we estimated it by assuming a coefficient of variance of 50% for OC and BC and of 80% for SO 2 . We assessed these coefficients of variance from all other observations providing standard deviation and mean. The samples for which the standard deviation was estimated are marked in Table S17.
-If the sample size for measurements was given as a range (e.g. "3 or 4"), we always took the lower number as N . When N was larger than 1 but not given, we assumed N = 2.
-Following Bond et al. (2004), we assumed that EF s follow a log-normal distribution. From the sample mean m and 5 standard deviation s, the mean µ and the standard deviation σ of the log-normal distribution were calculated: -For each emission sector, we randomly drew samples from the log-normal distributions with the calculated µ and σ.
-We used three different methods how to weight the different samples: i) every measurement (in Table S17) had the same 10 weight; ii) the measurements were weighted with N ; iii) every study had the same weight (differentiated by horizontal lines in Table S17). The three weighting methods can result in very different estimates, e.g. for OC emission factors from wood combustion (medians of 2.4 g kg −1 , 0.90 g kg −1 , and 3.0 g kg −1 , respectively).
-From the randomly generated samples, we calculated the median and the lower and upper quartile for each weighting method (the median is somewhat smaller than the expected value; Bond et al., 2004). The medians of the three weighting 15 methods were then averaged, and the same was done for the quartiles.
We considered the median to be the intermediate estimate for EF and the quartiles to be the lower and the upper estimates.
We did not choose more extreme percentiles for the lower and the upper estimates because the large variability in the measurements of EF reflects the high variability in burning conditions (e.g. smoldering versus flaming), fuels, combustion devices, and measurement devices. We explicitly wanted to consider these different conditions and not only sample from one (or a few) 20 measurements conducted under specific conditions. In general, more measurements for BC and OC than for SO 2 are available; however, since SO 2 emissions from biomass burning are small, we do not consider this to be an issue. In Sect. S3, the estimated EF s for the different sectors as well as the weighting of these different sectors are discussed. As described there, we estimate that 20% of the fuel consisted of agricultural waste, 40% of charcoal (in terms of wood that needs to be converted to charcoal), and 40% of wood. The combined aerosol emission factors were thus calculated as: For the intermediate scenario, we inserted the medians for EF agr , EF chw , and EF wood . For the low and the high estimates, the lower and the upper quartiles were used, respectively. The values for EF combined can be found in Table 1.

Aerosol emissions from crop residue burning
In the Greco-Roman world, fire was widely employed to fertilise fields, to create or regenerate pastures, to control pests, or to hunt (Ascoli and Bovio, 2013). In all simulations where we included the impact of anthropogenic aerosols, we considered the impact of humans on fires by (i) reducing the (natural) fire emissions calculated from CBALONE-SPITFIRE to account for crop and pasture areas (Sect. 2.4) and (ii) estimating fire emissions from pasture and crop residue burning. In the following, we 5 will first estimate aerosol emissions from crop residue burning before deriving estimates for pasture burning.
Next to fallow, crop rotation, and green manuring (i.e. plowing legumes in the soil; White, 1970), burning crop residues on the field was one method to increase the fertility of soils in Roman agriculture (Spurr, 1986). Since crop residues are burnt after harvest, emissions from open crop residue burning have a strong seasonal cycle (as for example shown for present-day China by Li et al., 2016b;Zhang et al., 2016). Today, harvest of cereal in the Mediterranean region takes place approximately 10 from the beginning of May to the end of August 2 , which is in accordance with summer being the time of harvest in Roman Italy (Spurr, 1986). We thus spread the fire emissions from crop residue burning over these four months; we neglected that a part of the harvest might have been burnt after August due to drying in the field. Since fires can get out of control at very high temperatures and are unlikely to be ignited at very low temperatures (Pfeiffer et al., 2013), we assume that crop burning cannot occur when the monthly surface temperature averaged over 20 years is below 0 • C or above 30 • C. The emitted mass 15 for these months without burning was shifted to the other months if there were any. The regions where temperature exceeded these thresholds were calculated offline using the surface temperature simulated with ECHAM-HAM-SALSA (without human impact; climatological analysis).
The aerosol emission fluxes [kg m −2 s −1 ] were calculated with the following equation (adapted from Webb et al., 2013): Cereals were the nutritional basis in the Roman Empire (Witcher, 2016) and typical yields were recorded by the agronomists of Classical Antiquity. Based on ancient sources and the work of Goodchild (2007) and Hopkins (2017), we estimated that 30 yields representative for the whole Roman Empire are in the range between 500 kg ha −1 and 1000 kg ha −1 (Sect. S4). was consumed per person per day. We roughly estimated that the wheat production was around 50% higher than the wheat consumption (resulting in 1.2 kg of crop yield per person per day) to consider that a part of the produced crop yield was lost 5 through, e.g. transport or insect damage (Spurr, 1986), used for fodder (Spurr, 1986), or needed for seeding (Hopkins, 1980).
For the low emissions scenario where the population size originates from HYDE11, this wheat production results in a yield of 440 kg ha −1 with the crop area based on HYDE11. For consistency, we use the simulated crop areas from ECHAM-HAM-SALSA (which are somewhat smaller; Sect. 2.3) for this calculation instead of the original input data. This yield is a bit lower than the estimates of crop production mentioned above (500-1000 kg ha −1 ). However, considering that part of the crop area 10 was used for other purposes than planting cereals or legumes (e.g. vineyards, which were not burnt) and that the crop area estimates from HYDE11 include areas that can lie fallow, the number seems reasonable: if we assume that fallow took place every second or third year (Goodchild, 2007), "actual" crop yields increase from 440 kg ha −1 to ≈ 660-880 kg ha −1 . For the intermediate scenario where population size is larger than in the low emission scenario, the HYDE database results in high but still reasonable estimates of crop yield (660 kg ha −1 ; ≈ 990-1320 kg ha −1 if fallow is considered). The high population 15 density in the high emission scenario would require an unrealistically high crop yield when combined with the crop area of HYDE11 (1100 kg ha −1 ; 1650-2200 kg ha −1 with fallow). When using the KK11 land use reconstruction, which has a much larger crop area, the needed crop yield would be 430 kg ha −1 (≈ 645-860 kg ha −1 ), which is regarded as realistic. Therefore, the HYDE11 land use is chosen for the low and the intermediate scenarios while KK11 is chosen for the high scenario.
For the yields, we took the values mentioned above (Table 2), which are calculated based on the crop area and the population size, for our calculations. Since our approach assumes a constant wheat consumption per person per day across all scenarios, we assure that crop yield provides in all cases the necessary food to feed the entire population. 5 In Sect. 2.5, we assumed that 20% of the fuel consisted of agricultural waste. Therefore, we should account for the fact that the higher the assumed population in our scenarios, the more crop residue is taken from the field. On the one hand, part of this agricultural waste used as fuel consisted not of cereal crop residue but e.g. of dung or olive pits, which speaks for a larger p b .
On the other hand, also other purposes of crop residue than fuel depend directly on the population density, e.g. the residue mass that was used for construction purposes or for filling mattresses (Spurr, 1986), which speaks for a lower p b . Here we simply 10 assume that the crop residue taken from the field per person is equal to 20% of the consumed fuel mass, and thus arrive at p b = 0.87, p b = 0.74, and p b = 0.56 for the low, the intermediate, and the high scenarios, respectively.

Fraction of crop burnt Fr crop_burnt
Part of the remaining residue on the field was burnt. We assume that the fraction of residue burnt in the field did not depend on the remaining residue mass per area but rather on cultural practices. There are several options for what happened after 15 harvest with the remaining residue on the field: when chaff was short, straw could be used as livestock feed (Spurr, 1986) or the residues left could be directly grazed by ruminants (Spurr, 1986), as it happens also nowadays e.g. in Mediterranean North Africa (Yevich and Logan, 2003). However, Spurr (1986) states that animals rarely eat stubble but rather the edible weeds and grass which grow among it. Farmers that had no further use for the stubble left after harvesting, burnt it (Spurr, 1986). In ancient sources, burning was mentioned as method to destroy weeds or to fertilise the soil (Spurr, 1986). Spurr (1986) assumed 20 that especially small farmers who could not store straw and had few animals to feed would burn the stubble.
Nowadays, crop residue burning is no longer widespread in developed countries, and many countries in western Europe even forbid open field burning (Yevich and Logan, 2003). To have a rough indication for how much of the remaining residue might have been burnt in the Roman Empire, we used present-day estimates from developing countries as an indication (including countries that cultivate rice despite the fact that rice was not grown in the Roman Empire). Yevich and Logan (2003) found 25 that in 1985 of the available residue, 1%, 23%, and 38% was burnt in the fields in China, in India, and in the rest of Asia, respectively. A more recent estimate for crop residue burning in China arrives at much higher fractions than the study by Yevich and Logan (2003): based on satellite data, Li et al. (2016a) (Yang et al., 2008). In line with this, Li et al.
(2016b) estimated that the fractions of crop residues burnt in fields increased from 5% in 1990 to 23% in 2013. Another reason why crop residue burning has increased might be that the use of biofuels has decreased in rural China (Yang et al., 2008).
We do not expect that mechanisation generally leads to enhanced crop residue burning; modern-day technology can also be used to prepare fields for the next crop planting after harvest (Pfeiffer et al., 2013), which makes the burning of crop residue unnecessary.
In the Philippines and Indonesia, up to 65% and 73% of the residue, respectively, was burnt in fields in 1985 (Yevich and Logan, 2003). This is similar to the satellite-derived fraction of rice fields under burning for the Punjab region in India in 2015 (66%, PRSC, 2015); as reasons for the large burning in the Punjab region, highly mechanised farming, poor storage facility 5 for the straw, and lack of market demand for further use are mentioned among others (PRSC, 2015). A survey in Bangladesh showed that only 3% of the farmers burn the total residue but that 37% burn the lower part of it (Haider, 2011); the surveyed farmers manually gathered the residue from the field.
By the use of p b , we have already accounted for the part of residue used as biofuel. If we remove this part of the residue from the estimates above, then the fractions of burning in the field would increase. Based on all the studies that we have mentioned,

Other variables in Equation 5
Modern values of the harvest index (= the ratio of grain yield to the total plant mass) for wheat range between 0.4 and 0.6 (Hay, 1995). The harvest index was typically lower (around 0.3) at the end of the nineteenth century (Hay, 1995;Sinclair, 1998), but the sometimes exceptionally high yields in ancient times could indicate that the harvest index might have been higher at 20 this time (Sinclair, 1998). For all simulations, we used a harvest index of 0.35, which corresponds to s = 1.9. Following Webb et al. (2013), we chose d = 0.85 and C f = 0.9. The emission factors for burning of crop residues on the field were estimated using the same method as described in Sect. 2.5.3 using values from Table S17 (key 3).

Aerosol emissions from pasture burning
According to the agronomist Columella (who lived in the first century AD), pasture from long fallow was burnt in late summer 25 to achieve more tender growth (Spurr, 1986 where Fr pasture is the fraction of the total gridbox covered by pasture, Fr pasture_burnt is the fraction of pasture that is burnt per time [s −1 ], F stands for fuel biomass consumption [kg dry_matter m −2 ; i.e. the amount of fuel that is actually burnt], the emission factors because these variables have large uncertainties.

Fuel biomass consumption F
The increase in European grassland productivity over the last decades has been small compared to crop (Smit et al., 2008). The spatial variability of grassland productivity is quite large within Europe, ranging from ≈ 0.15 kg m −2 yr −1 in the Mediterranean region up to 0.65 kg m −2 yr −1 in the Atlantic zones; the median over the different climate zones of Europe is 0.33 kg m −2 yr −1 10 (Smit et al., 2008). As a consequence, we expect that also the fuel load and the fuel biomass consumption show spatial variability. In contrast to fuel biomass consumption, fuel load refers to fuel that is available (but not necessarily burnt) and often refers to aboveground biomass only. Nevertheless, values for fuel biomass consumption and fuel load are quite comparable because grass easily burns and mainly aboveground biomass is consumed during pasture burning.
Values of fuel loads (kg dry_matter per area) frequently lie in the range between 0.3 and 0.5 kg m −2 for pastures and grass-  Based on these studies, we estimate that a reasonable average number for fuel biomass consumption is F = 0.35 kg dry_matter m −2 .

Fraction of pasture burnt Fr pasture_burnt
In the Mediterranean region, pasture burning has continued throughout history to the present day and is one component that has shaped the diversity of Mediterranean landscapes as we know them (Montiel and Kraus, 2010). In temperate and boreal Europe, the burning of grasslands for pasture was common on lands which were too poor in nutrients for agriculture, e.g. 5 heathlands (Montiel and Kraus, 2010).
In general, the abundance of prescribed burning depends on the accumulation of biomass: the higher the accumulation, the shorter the fire interval is. As a consequence, the fire interval depends on rainfall and grazing pressure (Weir et al., 2013;Frost and Robertson, 1987), thus showing pronounced regional variability. In the following we summarise guidelines for the rate of prescribed burning from different regions around the world. 10 For phryganic rangelands in Greece, it is recommended to set fire every 3 to 4 years, which allows to have a good herbage production and at the same time to suppress undesirable dwarf shrub (Papanastasis, 1980). In South Africa, Oluwole et al. (2008) found that the recovery period should be 3 years for optimum productivity in the absence of grazing. In line with this, found that grass richness, evenness, and diversity was high for sites with high rainfall when frequent burning was applied in the dry season (1-to 3-year return intervals), whereas Little et al. (2015) conclude that annual burning combined with intensive grazing has a detrimental effect on plant species diversity and structure. In Australia, single fires caused a short-term reduction of yield and cover of pastures in the following year, but fast recovery occurred for most burning regimes. However, perennial grasses were reduced at the expense of annual grasses, which is why burning every 5-6 and 4-6 years for arid short grass and 20 ribbon blue grass, respectively, are recommended (Dyer, 2011). This is in agreement with the findings of Norman (1963Norman ( , 1969 for native pasture on Tippera clay loam in the Katherine region. For North America, the recommended fire-return-interval of prescribed patch burning is 3 years in areas with rainfall above ≈ 760 mm per year and 4 years in drier regions (Weir et al., 2013).
One could argue that farmers in the past did not necessarily follow these present-day guidelines. However, traditional 25 knowledge of prescribed burning has been lost in many European areas (Montiel and Kraus, 2010). Guidelines thus partially re-establish knowledge that our ancestors had. On the one hand, if burning every year reduces the productivity of many grasslands, we think that it is unlikely that ancient farmers burnt their fields annually. On the other hand, a too long period without burning is also unlikely since this e.g. allows the growth of unwanted species and can have adverse effects on the ecosystem (SANBI, 2014;Papanastasis, 1980). According to the summarised literature, burning every ≈ 3 years seems to be 30 a reasonable intermediate estimate. Therefore, we assumed that 30 % of the pasture area is burnt per year for the intermediate scenario.
For the low and the high emission scenarios, we changed this fraction by a factor of 2 and thus arrive at 15 % and residue -that no pasture was burnt in months which are very cold or hot (monthly average temperatures below 0 • C or above

Emission factor EF pasture
Using the method described in Sect. 2.5.3 and the measurements from Table S17 (key 4), the emission factors in Table 3 were derived. We again used the lower quartiles for the low, the medians for intermediate, and the upper quartiles for high emission 5 scenarios.

Simulations
We conducted six time-slice experiments representative for the period around AD 100 (Table 4) In the standard ECHAM-HAM-SALSA version, the cloud droplet number concentration in a cloud cannot be lower than 40 cm −3 (or optionally 10 cm −3 ; Tegen et al., 2019). For all simulations of this study, this minimum CDNC was lowered from 40 cm −3 to 1 cm −3 (and the model was retuned) because the aerosol concentrations -and therefore most likely CDNCs -were considerably lower around AD 100 than today. 25 3 Results

Climate impact of anthropogenic land cover change
The simulations LCC_HYDE and LCC_KK were compared with the simulation no_human to quantify the impact of anthropogenic land cover change. Tables S8, S9 show seasonal and annual mean changes in some climate variables, including all that are discussed below.
30 Table 4. Overview of the different simulations. Note that the anthropogenic land cover change is applied in ECHAM-HAM-SALSA, but not in CBALONE-SPITFIRE, which was used to calculate aerosol emissions from natural fires exclusively. The replacement of natural vegetation by crop and pasture leads to a decrease in SOA precursors, which is most likely responsible for the significant 3 changes in CDNC, LWP, and CRE for LCC_KK (Table S9; all cloud properties are grid means).
In our simulations, the surface albedo can either increase or decrease significantly depending on the season and the region.
Deforestation can introduce brighter vegetation types and thus an increase in surface albedo (e.g. in Spain; Fig. 2a,b). However, especially after harvest, crop has a smaller canopy area fraction than natural vegetation. Thus, more of the soil, which is in 5 some regions darker than the natural vegetation, is exposed to radiation. As a consequence, the annually averaged surface albedo decreases in parts of Europe (Fig. 2a,b) because the surface albedo there increases in summer but decreases in other seasons (not shown).
Land cover change also alters latent heat (LH) and sensible heat (SH) fluxes. The annual mean evaporative fraction (Evap_frac= LH SH+LH ; only calculated if both SH and LH are upward fluxes, otherwise set to zero) shows regionally no significant changes 10 ( Fig. 2c,d). In contrast, regional changes in the annual mean turbulent flux (F turb = SH + LH; again only calculated if both SH and LH are upward fluxes) are significant for KK11 (Fig. 2f).
For HYDE11, the land surface temperature shows hardly any significant changes (Fig. 2g). For KK11, our results suggest that the changes in turbulent flux have a larger influence on the land surface temperature than the changes in surface albedo and evaporative fraction (Fig. 2). The decrease in turbulent flux regionally leads to a less efficient heat transport from the surface to 15 the atmosphere, thus enhancing the land surface temperature (Fig. 2h). In contrast, the land use in Smith et al. (2016) induces a distinct cooling over Europe around AD 1. We think that the opposite signal is mainly caused by differences in the vegetation or atmospheric models; intercomparison studies have shown that not all models agree on the sign of annual mean temperature changes induced by deforestation over North America and in mid-latitudes (Lejeune et al., 2017;Winckler et al., 2018). In our simulations, significant warming occurs south of 40°N, which is in qualitative agreement with present-day satellite studies (Li 20 3.2 Anthropogenic aerosol emissions around AD 100 and their magnitude compared to natural fire emissions For the majority of scenarios and aerosol (precursor) species, the emissions from pasture burning contribute most to the total anthropogenic emissions (Fig. 3). The emissions from fuel consumption are smaller in most cases but on the same order of magnitude (Fig. 3). Averaged over the whole year, the emissions from crop residue burning are relatively small; however, in summer, the emissions are in some cases comparable to pasture burning and/or fuel consumption (e.g. Fig. 3b).   (a)   Anthropogenic ACCMIP emissions in AD 1850 include the sectors industry, land transport, maritime transport, residential and 5 commercial combustion, agricultural waste burning on fields, and waste. Not all of these sectors produced aerosol emissions in AD 100; our emissions thus only include fuel consumption (both due to industrial as well as residential combustion), agricultural waste burning on fields, and pasture burning. Emissions from pasture burning are not an own ACCMIP sector; we expect that the anthropogenic ACCMIP emissions in AD 1850 could be somewhat higher if pasture burning were included.
The BC emissions from the high scenario in AD 100 are nearly as high as the emissions in AD 1850, whereas the emissions  (a) (b) (c) Figure 6. The same as Figure 5 but showing per capita emissions. in some cloud properties and radiative effects, averaged over the study domain and the whole year, are shown in Table 5.

5
Tables S10, S11, and S12 include more variables and show also seasonal averages for the low, the intermediate, and the high emission scenarios, respectively.
The BC and OC burdens (= vertically integrated mass) change significantly for all emission scenarios (Tables S10 to S12), with stronger changes for BC, whereas the SO 4 burden does not change significantly. This is due to the different emission sectors of the three aerosol types: while fire emissions are the only natural source of BC particles, OM also has biogenic 10 sources, and SO 4 predominantly forms from gas-to-particle conversion of volcanic and oceanic precursor emissions. The changes in aerosol burdens show large seasonal differences: in summer, the BC and OC burdens decrease for the low and the intermediate emission scenarios because of the reduction in natural fires (Sects. 2.4,3.2). In winter, the BC and OC burdens increase pronouncedly for all emission scenarios due to anthropogenic aerosol emissions.
Averaged over the study domain and the year, ARE changes by 0.06 W m −2 , 0.12 W m −2 , and 0.12 W m −2 for the low, 15 the intermediate, and the high emission scenarios, respectively (Table 5). Largest increases occur over North Africa and Europe ( Fig. 7a-c). In contrast to ARE, the changes in CRE are negative and considerably more pronounced: −2.11 W m −2 , −4.17 W m −2 , and −7.63 W m −2 for the three scenarios, respectively. The changes predominantly occur over Europe (Fig. 7df). The cooling effect of clouds is enhanced because they become optically thicker and more abundant due to aerosol-induced increases in CDNC and LWP (Tables S10 to S12). Overall, the negative changes in CRE are larger than the positive changes 20 in ARE. As a consequence, the anthropogenic aerosol emissions induce a decrease in land surface temperature (Fig. 7g-i).
Aerosol particles are also relevant because of their impact on human health. According to Makra (2015), air pollution had substantial consequences in cities in ancient times. As an example, the philosopher Seneca wrote that he needed to leave Rome in order to escape the smoke and the kitchen smells and to feel better (Makra, 2015). As maintained by the World Health Organisation (WHO;, the annual mean PM2.5 concentration should not exceed 10 µg m −3 . Figure 8 shows the background surface PM2.5 concentrations as well as the increases due to the anthropogenic emissions. For the low emission scenario, changes in the PM2.5 concentration are below 1 µg m −3 everywhere. For the intermediate emission scenario, significant increases mainly range between 0.1 and 2 µg m −3 , which is for most places insufficient to increase the background concentrations 5 to values above 10 µg m −3 . Only for the high emission scenario the increase in PM2.5 concentration is pronounced enough to exceed the WHO threshold in larger areas (e.g. Algeria, Egypt). The maximum PM2.5 concentrations are to some degree underestimated due to our coarse model resolution. As an example, Li et al. (2016c) found that the maximum simulated PM2.5 concentrations decrease by -21% if the resolution decreases from 0.5°× 0.66°to 2°× 2.5°.
(a)  Because ECHAM-HAM-SALSA is an atmosphere-only model, SST and SIC were prescribed. Therefore, temperature does not react to a perturbation as much as with an interactive ocean, implying that, in the absence of attenuating feedbacks, our simulated changes in temperature should be underestimated. We repeated two of the simulations (no_human, LCC_HYDE_int) with a mixed-layer ocean (MLO) model having a depth of 50 m (≈ 50 years of simulation, half of it model spin-up), and arrived 5 indeed at larger changes in land surface temperature (4 times stronger decrease in land surface temperature averaged over our study domain). Nevertheless, the temperature patterns occurring with fixed SST remain visible with the MLO setup: the land surface temperature decreases mainly over Europe due to anthropogenic aerosol emissions (Fig. S4). As expected, the aerosolinduced cooling is more widespread with the MLO, but it is basically limited to the Roman Empire. Since the changes in variables affecting temperature (e.g. CRE) are comparable for the different model setups, qualitatively similar, though stronger 10 changes in the land surface temperature occur with the MLO. Nevertheless, we admit that a more detailed analysis using the MLO, including all scenarios, will provide more quantitative and perhaps partly different results. However, this is beyond the scope of this study.
The two land cover reconstructions that we chose differ widely and thus provide a measure for uncertainty. However, also the soil albedo and vegetation cover of the simulation no_human, which have a strong influence on the climate response, are 15 subject to uncertainty. We found that the contribution of grasslands to the natural vegetation is quite high in our simulations. As an example, grasslands and tundra contribute ≈ 45% to the natural vegetation between 0°E and 30°E and 30°N and 60°N.
If the forest fraction were underestimated in our simulation, the impact of the anthropogenic land cover change would likely also be underestimated. This would especially be the case for creating pasture, since grasslands are preferentially converted to pasture in the model (Reick et al., 2013). Furthermore, the crop phenology of JSBACH, which considers both winter and 20 summer crop, allows two harvests of summer crops during one year in the extratropics depending on the heat sum. Generally, the productivity of crops was smaller in the Roman Empire than today, with fallow every two to three years. The number of harvests is thus likely overestimated in JSBACH (but not in the offline calculated crop residue burning emissions), which affects for example changes in surface albedo (which could be more positive in reality).
The anthropogenic aerosol emissions that we calculated are only rough estimates. To calculate aerosol emission factors for 25 fuel consumption, we used measurements from present-day fireplaces and traditional stoves. The measurements show a large variability, which is caused e.g. by the fuel moisture, the burning device, the fire conditions (flaming versus smoldering), or the instrumental setup. To account for the large range of observed emission factors, we considered many measurements. Nevertheless, we cannot rule out that typical stoves used in the Roman Empire had systematically different emission factors than the majority of burning devices that we considered. Furthermore, the uncertainties of important parameters such as population 30 size are considerable when going so far back in time. The emissions are thus highly uncertain, and it is extremely challenging to quantitatively verify them with palaeo records. At least there is indirect evidence for aerosol emissions associated with fuel consumption and agricultural burning: marble turned grey in antique towns due to smoke, and laws were introduced against air pollution (Makra, 2015). Moreover, pollen and charcoal records show high positive correlations between fire and crops, weeds, and shrubs in Mediterranean and temperate regions in and around the Alps between 2300 BC and 800 AD (Tinner et al., 2005(Tinner et al., , 2009. Palaeo records further suggest that controlled burning was used to introduce and establish sweet chestnut in some regions (Morales-Molino et al., 2015).
Another simplification is our lack of spatial and temporal variations: for many variables, we estimated one "typical" value for the whole Roman Empire, and our anthropogenic emissions show no trends over the years (e.g. caused by wars). The time 5 around AD 100 was a relatively stable period, characterised by expansion of infrastructure, economic wealth, and quite low military activities. Nevertheless, the emissions were more dynamic in reality than in our simulations, e.g. due to Trajan's Dacian Wars (AD 101/102 and AD 105/106), which caused population movements as well as a possible change in anthropogenic aerosol emissions (e.g. regional emissions associated with warlike activities such as burning of villages).
The impact of our simulated anthropogenic aerosol emissions strongly depends on the natural background and its season- be overestimated. On the one hand, this could indiciate that our background aerosol concentration is potentially overestimated as well. On the other hand, our model neglects pure biogenic aerosol nucleation and nitrate aerosols. Using the GLOMAP model, Gordon et al. (2016) have shown that the global annual mean CCN concentration (at 0.2% supersaturation and at cloud base level) was 12% higher in the pre-industrial atmosphere when pure biogenic nucleation is considered. The effect is not the same over the year since terpene emissions are larger in summer (Gordon et al., 2016). Therefore, our simulated background 20 concentration when we see the highest impact of anthropogenic aerosols (winter, spring) might not be affected too much by the neglect of biogenic nucleation.
The aerosol background does not only depend on emissions, but also on other parameters such as the simulated size distribution or the calculation of removal processes. A simplification of our study is that most of the tropospheric chemistry is not calculated, hence we prescribe the oxidants that are important for the sulfur chemistry and the SOA formation. Also the ver-25 tical distribution of the aerosol particles is essential when their effect on radiation and clouds is analysed. In ECHAM-HAM, vertical mixing is strong in the lower troposphere (Veira et al., 2015), and thus the anthropogenic aerosol particles emitted near the surface can easily reach altitudes where clouds form. If this vertical mixing were overestimated in the model, then the aerosol-cloud interactions would likely be overestimated as well.
Last but not least, aerosol-radiation and aerosol-cloud interactions from aerosol-climate models are uncertain. As an exam-30 ple, some studies show that the models typically overestimate the effect of aerosols on the cloud liquid water content, at least in some regions (Bender et al., 2018). The aerosol effective radiative forcing also depends on the minimum CDNC value (Hoose et al., 2009); since we lowered the minimum CDNC to 1 cm −3 , the aerosols have a stronger impact on radiation than with the standard setup of ECHAM-HAM-SALSA where the minimum CDNC is 40 cm −3 . We conducted additional simulations showing that the total ERF due to aerosols between AD 1850 and AD 2000 is ≈ −2.4 W m −2 with our decreased minimum 35 CDNC value. Our ERF ari+aci thus lies at the lower end of model estimates and is outside the range of IPCC's expert judgment (−1.9 to −0.1 W m −2 ; Boucher et al., 2013). Nevertheless, we expect that the anthropogenic aerosol emissions around AD 100 would still increase the CDNC in the Roman Empire (and thus induce a cooling effect) with the standard minimum CDNC value of 40 cm −3 : the simulated CDNCs are above 40 cm −3 in the lower troposphere for all seasons when averaged over our study domain (not shown). The strong total ERF ari+aci with a minimum CDNC of 1 cm −3 therefore mainly results from other 5 regions, e.g. the Arctic and the Antarctic.

Conclusions
The hypothesis of this study was that anthropogenic activities associated with land cover and aerosols already had a noticeable influence on some climate variables in the Roman Empire around AD 100. To test this hypothesis, we first created three different scenarios of anthropogenic aerosol emissions for the Roman Empire. These scenarios were combined with two existing land 10 use datasets.
Simulations with an aerosol-enabled global climate model showed that the anthropogenic land cover reconstruction by KK11 (Kaplan et al., 2011 induces significant decreases in turbulent flux and thus a warming in parts of North Africa and the Middle East, whereas the reconstruction based on HYDE (Klein Goldewijk et al., 2010 has nearly no impact. Mainly over Central and Eastern Europe, anthropogenic aerosol emissions lead to an enhanced cooling effect of clouds (i.e. a more 15 negative cloud radiative effect), the extent and strength of which depends on the different scenarios. Our model likely overestimates changes in the cloud radiative effect to some degree due to the lowered minimum CDNC (Hoose et al., 2009). Overall, land use thus generally increases the land surface temperature in our simulations, whereas anthropogenic aerosol emissions decrease it. Since we prescribe the SST, our simulated changes in temperature are strongly underestimated. Simulations with an interactive ocean are therefore needed for a more quantitative analysis. 20 Around AD 100 temperatures in Europe (as well as China and the Northern Hemisphere in general) were warmer than the average over the last two millenia (PAGES 2k Consortium, 2013;Ge et al., 2013;Christiansen and Ljungqvist, 2012), an era that has been called the "Roman Warm Period". While our results imply that anthropogenic land cover change may have regionally contributed to this warming, aerosol-cloud interactions would have attenuated it, suggesting other causes of the Roman Warm Period, e.g. ocean dynamics or solar forcing. 25 Our scenarios show that pasture burning could have been an important source of aerosol particles. A better understanding of the processes that drive the frequency, the seasonality, and the emissions of pasture burning could therefore be essential to quantify the anthropogenic impact far back in time. These processes could be implemented in a vegetation-fire model that is coupled with an aerosol model. Simulations with a higher spatial resolution would moreover allow to account for the large regional variations in aerosol emissions within the Roman Empire -information which can be provided by archaeologist and 30 historians. Therefore, further work should include collaborations with them in order to incorporate better understanding of human behaviour in Classical Antiquity, particularly on fuel consumption and the timing and extent of the use of fire. Such collaborations could also allow to assess the effect of economic crises (e.g. in the third century AD) and wars on aerosol emissions.
Code availability. TEXT

Data availability. TEXT
Code and data availability. TEXT