A virtual water network of the Roman world

The Romans were perhaps the most impressive exponents of water resource management in preindustrial times with irrigation and virtual water trade facilitating unprecedented urbanisation and socioeconomic stability for hundreds of years in a region of highly variable climate. To understand Roman water resource management in response to urbanisation and climate variability, a Virtual Water Network of the Roman World was developed. Using this network we find that irrigation and virtual water trade increased Roman resilience to interannual climate variability. However, urbanisation arising from virtual water trade likely pushed the Empire closer to the boundary of its water 20 resources, led to an increase in import costs, and eroded its resilience to climate variability in the long-term. In addition to improving our understanding of Roman water resource management, our cost-distance based analysis illuminates how increases in import costs arising from climatic and population pressures are likely to be distributed in the future global virtual water network.


Introduction
Trade is central to safeguarding food security under the twin pressures of growing demand and intensified climate variability (Godfray et al., 2010;Schmidhuber and Tubiello, 2007).The redistribution of food through trade sustains populations where local food resources are insufficient to meet demand or where climatic variability causes low yields ( Schmidhuber and Tubiello, 2007).Trade in food is intimately linked to the freshwater resources of trading regions with up to 90% of human freshwater use going to agricultural production (Hoekstra and Chapagain, 2008;Shiklomanov, 2000).The freshwater resources embodied in food production and traded among regions is known as virtual water (VW) (Allan, 1998) and by tracking VW flows it is possible to quantify how freshwater resources are redistributed around the globe (Hoekstra and Chapagain, 2008).Great strides have been made to empirically describe the global trade in VW (Carr et al., 2012a;Dalin et al., 2012a;Konar et al., 2011;Suweis et al., 2011) and quantify the volume of VW flows among regions (Hanasaki et al., 2010;Shi et al., 2014;Suweis et al., 2013).Studies have shown that VW predominantly flows from regions with a surplus in water resources (water rich) to those with insufficient resources to meet local demand (water poor) ( Dalin et al., 2012a;Liu and Savenije, 2008;Shi et al., 2014).However, Konar et al. (2011) found that while VW redistribution saves water on average globally, many bilateral VW trade links are irrational from a water savings perspective and exist instead for complex socioeconomic reasons such as trade agreements, wealth disparity, agricultural subsidies and so on (de Fraiture et al., 2004).As a result, isolating the impact of climate variability and population demand on VW flows is challenging because the imprint of complex socioeconomic forcings overprint and are intertwined with the response of the VW network to climate and population forcings (Dalin et al, 2012b;Suweis et al., 2011).Additionally, the complexity of socioeconomic forcings and crop response to climate change make future predictions on VW trade and VW content highly uncertain (Fader et al., 2010;Konar et al., 2013;Zhao et al., 2014).Sivapalan et al. (2012) recommend studying past society's relations with water, a term they refer to as historical socio-hydrology, to understand fundamental processes linking humans and water resources.They propose that water has played a role in the growth, evolution and eventual collapse of many past societies.Therefore by studying the relation of past societies with water we can illuminate factors that are also relevant in present-day, such as how close we are to reaching the planetary boundaries of current fresh water resources (Bogardi et al., 2013;Rockström et al., 2009).Historical socio-hydrology principals were applied in an analysis of the coevolution of society and water resources over the last 2,000 years in the Tarim river basin of Northwest China (Liu et al., 2014).Based on historical reconstructions, Liu et al. (2014) developed a conceptual model which identified a positive, destructive feedback loop related to economic gains from water resources and a negative, restorative feedback loop related to the impacts of exploitation of water resources.That same model was shown to be representative of the changes in the Murrumbidgee catchment in Australia over the last 100 years indicating that an historical basis for understanding socio-hydrological systems shows promise (van Emmerik et al., 2014).In a modelling context, one of the principal advantages of using historical data is that it is possible to constrain simulations based on historical reconstructions and compare model output with what actually happened (given the inherent uncertainties of historical reconstructions) (Cornell et al., 2010;van der Leeuw et al., 2011).Such an approach can indicate if the assumptions and processes incorporated in models are valid.A modelling approach can also be revealing in an historical context (Cornell et al., 2010).For example, the use of physically-based models in an historical context also facilitates the refinement of theories about past human relations with their environment based on what was physically possible given the constraints of the environment.
In the preindustrial period, the Roman Empire were likely the greatest exponents of virtual water trade as evidenced by the widespread trade in water resources, particularly grain, throughout the Mediterranean and Black Sea region (Erdkamp, 2005;Kessler and Temin, 2007;Rickman, 1980;Scheidel, 2010).Supplying the main cities of the Empire with sufficient grain was one of the principal preoccupations of the ruling elite throughout the lifetime of the Republic and Empire, to the extent that a stable supply of grain to the city of Rome became personified by the deity Anonna (Mazoyer and Roudart, 2006;Rickman, 1980).In a close parallel to current demographic trends (Chen, 2007;United Nations, 2012), an explosion in urban populations during the Late Republican era (Bowman and Wilson, 2011) led many cities to overshoot their local ecohydrological carrying capacities bringing about an increased reliance on imports of VW (Erdkamp, 2005).Similar trends, but of a much greater magnitude, are seen present day across the globe in countries such as China where rapid urbanisation, increased affluence and relaxing of trade restrictions have brought about a 20 fold increase in VW imports in less than a quarter of a century (Hubacek et al., 2009;Shi et al., 2014).
As with current society, the Romans sought to secure food security in two principal ways: through in situ water resource management using rainfed agriculture and irrigation (Fader et al., 2009;Torell et al., 1990;Wada et al., 2011) and through VW trade (Barnaby, 2009;Shi et al., 2014;Yang and Zehnder, 2001).Irrigation enabled the Romans to maximise exploitation of local water resources whilst VW trade allowed them to inhabit regions where local water resources were insufficient for the resident population (Garnsey, 1998;D'Odorico et al., 2010).The Romans also made use of large municipal grain stores which were replenished after each harvest owing to spoilage.These municipal stores acted as a buffer for when imports became disrupted (Erdkamp, 2005).Temporal market speculation on grain through hoarding is thought to have been limited in the Roman Period, however.Market speculation was a high risk venture owing to the loss in value of grain as a result of storage and the high uncertainty 90 associated with predicting surpluses or deficits in subsequent years (Erdkamp, 2005).As a result VW distribution predominantly responded directly to yield surpluses and deficits integrated over a short number of years rather than complex economic dynamics arising from speculation (Erdkamp, 2005;Horden and Purcell, 2000).
In terms of in situ water resource management the Romans made use of a wealth irrigation technologies such as dams, aqueducts, canals, cisterns, water wheels and Qanats (Barker, 1996;Wilson, 1997).The maintenance and operation of irrigation infrastructure was tightly controlled with users taxed on the extent of land they irrigated in regions such as Egypt or on the magnitude of their harvest in Spain, Sicily and Sardinia (Beltrán Lloris, 2006;Erdkamp, 2005).The Romans were far from the first Mediterranean civilisation to use such water management technologies but the extent and organisation was unprecedented and enabled them to achieve high agricultural yields 100 throughout their Empire (Barker, 1996).Not only did irrigation increase grain production but surface water was a far more reliable source of water compared with precipitation, particularly in large river basins such as the Nile delta, the Po Valley and the Orontes, Ebro and Vera catchments in present-day Syria and Spain (Beltrán Lloris, 2006;Butzer et al., 1985;van der Leeuw, 1998).
VW redistribution during the Roman period was comparatively simple compared with present day global trade in water resources (Erdkamp, 2005;Konar et al., 2011).Within the Roman Empire few artificial trade barriers existed, instead the redistribution of VW was driven by satisfying demand of urban centres from regions with a surplus principally by means of tributary donations (tax-in-kind in the form of grain) (Erdkamp, 2005;Scheidel, 2010;Temin, 2012).There was a parallel free market, however it is thought the urban grain supply was too important to be risked on the free market and was ensured by hierarchical methods (Rickman, 1980).The principal constraint on VW redistribution in the Roman Period was the 'struggle against distance' (Braudel, 1995).However, advanced shipping technology during the Roman Period combined with the relative safety of summer maritime travel within the Mediterranean facilitated unprecedented trade in bulk goods such as grain (Houston, 1988).As with present-day, trade costs of these bulk goods co-varied with distance (Hummels, 2007).However, transport by ship was significantly cheaper compared with overland transport owing to the difficultly in land-based transport of bulk goods by horse and cart (Jones, 1986;Scheidel, 2013); a feature of trade in bulk goods that remains despite modern advancements in transport technology (Limão and Venables, 2001).
In this paper we set out to examine how irrigation and VW trade contributed to Roman resilience against the twin pressures of urbanisation and climate variability.In order to understand this we have developed a Virtual Water Network of the Roman World.Our VW network contains two principal components: a hydrological model and a dynamic, agent-based VW redistribution network.We simulate yields under variable climate conditions using the hydrological model PC Raster Global Water Balance Model (PCR-GLOBWB) ( van Beek and Bierkens, 2009;van Beek et al., 2011).VW trade is simulated using Orbis, the Stanford Geospatial Network of the Roman World (Scheidel, 2013) as our network structure, with link weights reflecting transport costs at 200 AD associated with the 'struggle against distance' (Braudel, 1995).We do not model the feedbacks between the trade dynamics and hydrology (similar to the socio-hydrological approach proposed by Liu et al., (2014) and Elshafei et al., (2014) for example); rather we apply a scenario-based approach using historical reconstructions and the physical hydrological model as constraints.Our analysis of the Roman water resource management not only adds to our understanding of that civilisation but also helps us to understand the fundamental processes underpinning VW trade in present-day (Sivapalan et al., 2012).

Methods
The schematic of our Virtual Water Network of the Roman World is shown in Fig. 1.To summarise our methodology; we calculated yields using the global hydrological model PCR-GLOBWB based on estimates of the extent of Roman cropland in 200 AD from the History Database of the Global Environment (HYDE) (Klein Goldewijk et al., 2011).Land with a potential for irrigation was assigned within HYDE cropland regions based on the MIRCA dataset of Portmann et al. (2008).Natural landcover was assigned based on the Olson classification (Olson, 1994a(Olson, , 1994b)).The yield response to climate variability was calculated in PCR-GLOBWB with climate prescribed using meteorological observations over the period 1949-2000(Ngo-Duc et al., 2005).VW surpluses and deficits were calculated with VW demand based on HYDE gridded population estimates.Yearly VW surpluses and deficits were abstracted to Orbis and the redistribution of VW from VW rich to VW poor regions of the Roman Empire was simulated.A detailed description of our methodology follows.

Simulating cereal yields under variable climate
We computed VW based on cereals at 5' resolution under rainfed and irrigated cultivation using the Global Hydrological model PCR-GLOBWB (see van Beek et al. (2011) and van Beek and Bierkens, (2009) for a detailed description of the model).PCR-GLOBWB is a spatially explicit 3 layer (2 soil layers and 1 groundwater reservoir) hydrological model that computes the vertical water balance for different land cover types under prescribed meteorological conditions and routes the specific runoff to obtain discharge fields.One of the outputs of the vertical water balance is the actual transpiration (so-called green water) which was used here to estimate grain yields (Doorenbos and A.H. Kassam, 1979).When soil moisture is limiting, yield may be maximized for a healthy and fertilized crop if the crop water requirements are met by irrigation (blue water) and the crop transpires at the potential rate sustained by the atmospheric demand (Allen et al., 1998).Following this principle, yield can be taken to be proportional to the water use efficiency multiplied by transpiration (Zwart et al., 2010).The crop water requirements equal the difference between potential and actual evapotranspiration for the cropped area and correspond to the irrigation water demand when divided by the irrigation efficiency that accounts for conveyance and application losses.Using these principles, irrigation water demand and the realized yield were evaluated on a monthly scale with consideration of climate variability.To this end, the potential and actual evapotranspiration rates of cropped areas when fed by rainfall only were used to compute the irrigation water demand (see Wada et al., 2011 for details) and to ascertain what proportion of the irrigation water demand can be satisfied with the available discharge.
PCR-GLOBWB requires meteorological and land cover data as input.As meteorological forcing we used the National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) corrected by the Climate Research Unit (CRU) climate reanalysis dataset over the period 1949-2000(Ngo-Duc et al., 2005), which downscales NCEP/NCAR data to a regular 1-degree global grid with a daily resolution.Using current reanalysis data for the Roman period was deemed acceptable as the reconstructed Roman climate optimum was estimated to be comparable with the mean Northern Hemisphere temperature between 1961-1990(Ljungqvist, 2010)).
In terms of precipitation, early modelling studies had suggested that greater forest cover in the Roman period maintained a wetter climate (Reale and Dirmeyer, 2000;Reale and Shukla, 2000).However, historical, archaeological and paleoclimatological evidence indicates that the mean background climate in the Mediterranean during the Roman period was broadly similar to present day although there were likely centennial-millennial shifts in synoptic climate systems which would have made certain regions relatively drier or wetter on average at different times during the Roman Period (Büntgen et al., 2011;Dermody et al., 2012).The impacts of longer-term shifts in the synoptic climate systems on Roman water resource management will be assessed in a follow-up paper.The CRU TS 2.1 dataset only specifies variables for the global land mass so to ensure global coverage, the original NCEP/NCAR values were inserted if no values were specified.From this dataset, daily precipitation totals and the average temperature were used directly as model input.The model also requires reference potential evapotranspiration as direct input which was computed using the Hamon method (Allen et al., 1998), which only requires temperature as meteorological input compared to more complex equations.Monthly climatology's of wind speed and relative humidity were used indirectly to estimate the crop factors (see below).
To partition precipitation (rainfall, snow) into interception and throughfall and to prescribe the crop-specific potential evapotranspiration, PCR-GLOBWB requires the interception capacity, ground cover and the crop coefficient for each land cover type.The natural land cover parameterization is based on the Global LandCover Characterisation (GLCC) at 30" with the Olson classification (Olson, 1994a and1994b) and the parameter set of Hagemann (1999).
Irrigated areas were inserted using the MIRCA dataset of Portmann et al. (2008); (see van Beek et al., 2011;Wada et al., 2011 for details).The fraction of each cell assigned as crop and pasture land was defined based on History Database of the Global Environment (HYDE) reconstructions for 200AD at 5' horizontal resolution (Klein Goldewijk et al., 2011) (Fig. S3).HYDE does not explicitly account for crop rotation and the issue of crop rotation in the Roman period remains controversial with some authors claiming that no rotation was practised in Roman times whereas others claim two-field rotation was practised with one half of fields laying fallow at any one time (Fox, 1986).Based on White (1970), we adopt an intermediate value of continual three field cropping with 2 years of a cereal crop and 1 year of fallow assigned as sparse grassland according to GLCC.In irrigated regions we employ multi-season cereal cropping based on the crop calendars from the MIRCA dataset (Portmann et al., 2008).
The land cover parameterization is derived from the 30" distribution of the GLCC (Olson, 1994a(Olson, , 1994b)).In order to incorporate the information on cultivated area for the Roman period from the HYDE dataset, having a spatial resolution of 5', the distribution of cultivated and pasture areas was reconstructed at the resolution of 30".Within each 5' cell, all 30" cells were ranked on suitability; using the GLCC classification at 30", areas were delineated to represent respectively the presently cultivated areas and those under pasture.Within these areas, each cell was assigned a decreasing suitability with increasing slope owing to the Roman's preference for low lying, gently sloping land for agriculture (van der Leeuw, 1998).Outside the presently exploited areas, suitability was ranked using the slope parallel cumulative distance from the boundaries of these areas outwards.Suitability was then scaled between the minimum and maximum values to yield a range between 0 and 1.This suitability was then used to iteratively select the most suitable cells until the desired area was met.Precedence was first given to cultivated area, followed by pasture.The remaining area was filled with the reconstructed natural vegetation from the GLCC dataset.The resulting mosaic at 30" was consecutively used to compute the effective values of the land cover parameterization per land use type at 5'.Any remaining cells were assigned as semi-natural land cover types that were extrapolated spatially on the basis of the Holdridge Life Zones (Leemans, 1990(Leemans, , 1992)).For the semi-natural vegetation, a subdivision between short and tall natural vegetation was made on the basis of forest fraction.
Cropland was subdivided proportionally into irrigated and rainfed land on the basis of the MIRCA dataset, giving, with pasture, a total of five land cover classes within each cell.Monthly characteristics were prescribed to account for seasonal growth changes in cereals and natural vegetation.For short-natural, tall-natural and pasture land cover types, these values were based on the original Olson classification and the corresponding parameterization of Hagemann et al. (1999).For the irrigated and rainfed cropland, the crop factors and calendars were taken directly from the MIRCA dataset (Portmann et al., 2008) for cereals under rainfed and irrigated conditions and weighted by area.Water use efficiency for all cereals was assigned the value for winter wheat and crop yield taken to be equivalent to 25% of the total above-ground biomass compared to the 35% used by Zwart et al. (2010) for presentday crops, in line with estimates from Roman and pre-agricultural revolution sources (Erdkamp, 2005;Goodchild, 2007).It is important to highlight that we only calculated yields based on cereal crops whereas large portions of land would have also been given over to viticulture, olives, market gardens etc. (Columella, 70AD;Erdkamp, 2005).HYDE population values were used to calculate VW water demand as well as the workforce available for harvest.In addition to water availability, labour availability constrains the area that can be cultivated.The labouring population was calculated based on the grid-based population estimates from the HYDE dataset (Klein Goldewijk et al., 2011).
We estimated a harvesting period of 1 month with an average harvest area per person per day of 0.2 ha which equates to 6 ha per person per year.We restricted harvesting to the able-bodied population aged between 12 and 55.
Based on demographic life tables from Roman Egypt, this equated to 55% of the population that were capable of helping with the harvest (Frier, 1982).HYDE population values were also used to calculate VW demand based on a consumption of 200 kg of grain per person per year (Erdkamp, 2005).For the 20 most populous cities in the empire, the grid-based population values of HYDE were corrected using Chandler's (1987) estimates of Roman urban population.For each cell we subtracted the population demand from the realized yield providing yearly maps of surplus and deficit VW (Fig. S4).

Simulating virtual water redistribution
Orbis, the Stanford Geospatial Network of the Roman World forms the basis for our VW redistribution network of the Roman World (Meeks, 2013;Scheidel, 2013).Orbis broadly reflects the transport network in the Roman Empire around 200 AD with all links confirmed as Roman era transport routes although we cannot be certain that all were active in 200 AD.Orbis should be taken in the spirit for which it was intended, which is to outline the dramatic contrasts between terrestrial, fluvial, and maritime transportation expenses and the patterns they imposed on the flow of goods within the Roman Empire (Scheidel, 2013).Orbis contains a database of 751 roman towns and cities that form the nodes within our network.These cities are linked by (1,371 x 2) directed edge segments that represent the cost to transport a kilogram of grain in denarii along Roman roads, rivers and over sea in each month of the year based on Diocletian's edict on Maximum Prices and physical cost distance calculations (Scheidel, 2013).The links between each node have a cost representing transport in each direction.For example, up-river transport is more costly compared with down-river transport.We used transport costs for the month of June because the majority of grain was transported during summer months when sea conditions were calm (Erdkamp, 2005;Horden and Purcell, 2000).
We collapsed nodes within 10km of each other into 1 node owing to the resolution of the underlying gridded data, leaving us with 649 nodes.To simplify calculations in our dynamic redistribution model, we converted the directed network of Orbis into an undirected network by taking the average cost of the directed links between nodes resulting in a total of 1,371 undirected links.As we are interested in Mediterranean climate variability, we restricted our analysis to the part of the network that extends from 10W to 45E and 25N to 46N, however all simulations were carried out for the entire Empire.In order to convert grid-based surplus and deficit data to the Orbis network structure we assigned city regions using a Theissen polygon operation between our city nodes (see Fig. S1 for city regions).All gridded data within these regions were summed and applied to the relevant city or town node.Therefore certain nodes in the network were either VW rich or VW poor based on the total grain yield -total grain demand within that city region.The VW surplus and deficits in each node changed each year based on changes in yield owing to climate variability.We represent VW water imports and exports in terms of per person VW demand rather than cubic metres of water to make our findings more accessible to non-specialists in agronomy and hydrology.
Our VW redistribution network operates as a dynamical agent-based network (Wilensky, 1999).In line with our understanding of the Roman grain economy (Erdkamp, 2005;Scheidel, 2010), our network is demand driven with each VW poor node (nodes with a VW deficit) individually demanding VW from linked VW rich nodes.That is to say, that our network simulates a hierarchical grain supply system whereby urban centres ensure supplies through taxation in kind, constrained by the struggle against distance (Braudel, 1995).Although, much of the evidence we base our trade rules pertain to the city of Rome, it is likely that other major cities used similar methods to ensure grain supplies.Similarly to D'Odorico et al. ( 2010), we do not simulate VW trade between VW rich nodes although this may have occurred.Since the links in our network are undirected, flow direction is dictated by the VW potential among VW rich and poor nodes.Thus, VW flow in our network responds directly to changes in yields arising from climate variability.VW redistribution is simulated over 52 years of climate variability (Ngo-Duc et al., 2005), with a year ending when demand at all deficit nodes has been satisfied or all surplus nodes are depleted.We quantify the stress on the system in terms of the cost to import VW with costs measured at all VW poor nodes.

Yield response to climate variability
The yearly average simulated yield for cereals per 5' cell is shown in Fig. 2a with the contribution to the total from rainfed (Fig. 2b) and irrigated (Fig. 2c) land shown separately.The yields in kg ha -1 are shown in Fig. S5, however since HYDE cropland fractions vary per cell, the yield per 5' cell give a clearer impression of spatial variability in total yield amount.Rather than reporting VW partitioned into its green and blue component sources we partitioned VW into VW derived from rainfed and irrigated land.Yields from rainfed land derive only from green water whereas yields from irrigated land incorporate blue water where there is a shortfall in green water to meet the evaporative demand (van Beek et al., 2011).Our simulations indicate that the most productive rainfed agricultural regions are located in present-day Spain, France, the Po valley, Western Turkey and the Fertile Crescent (present-day Syria, Iraq and Israel) (Fig. 2b).Irrigation agriculture is also widespread (Fig. 2c), with the largest areas of irrigated agriculture located in Egypt, the Po valley, south-eastern Turkey, the Fertile Crescent and Spain.Rainfed agriculture accounts for 71.5% of the total yield in the region with irrigation accounting for the remaining 28.5%.The kg ha -1 yields (Fig. S5) are consistent with yield estimates based on Roman sources and yields prior to the agricultural revolution in Europe (Erdkamp, 2005;Goodchild, 2007).
Lower than expected yields are calculated for Sicily and present-day Algeria and Tunisia related to what is known from historical sources about the productivity of these regions (Erdkamp, 2005).The low yields in these regions are due to a probable underestimation of cropland fractions in the HYDE dataset (Fig. S3b).HYDE provides estimates of cropland fractions and population concentration at 5' spatial resolution globally for the entire Holocene using land suitability algorithms and back-calculating from current population and cropland distributions (Klein Goldewijk et al., 2011).Thus, it is not surprising that for certain regions cropland fractions are inconsistent with historical accounts for the specific date of 200 AD (Fig. S3) (Klein Goldewijk and Verburg, 2013).For the purposes of this paper it was decided to use unadjusted HYDE grid-based estimates of cropland to transparently show our methodology.
Proxy reconstructions indicate anomalously warm climate conditions during the Roman period owing to warm temperatures of the North Atlantic Ocean at the time (Bond et al., 2001;Desprat et al., 2003;McDermott et al., 2001).Fig. 3 shows the correlation between grain yield with temperature and precipitation averaged over land cells in the Mediterranean region (25N 46N and 10W 45E).Each point represents a single year of climate forcing with temperature and precipitation plotted on the x-axis and yield on the y-axis.Under warmer temperatures, grain yield significantly increased in both rainfed (p=0.001) and irrigated (p=0.001)regions (Fig. 3 a, c, e).Somewhat counterintuitively, yield significantly decreased in rainfed regions under increased precipitation (p=0.008) (Fig. 3b).

310
No significant relation was found between precipitation and yields in irrigated regions (p=0.62).Yield is calculated based on evapotranspiration, with warmer conditions bringing about higher evapotranspiration and thus higher yields where water is not limiting (van Beek et al., 2011).Yield decreases under increased precipitation owing to the negative relation between temperature and precipitation in most of the Mediterranean for the predominantly winterspring growing season (Fig. S6) (Portmann, 2008).Additionally, depending on soil type and average rainfall, transpiration can be limited in PCR-GLOBWB by oxygen stress in the soil caused by water logging ( van Beek and Bierkens, 2009).In irrigated regions there is no relation with precipitation because much of the growing period in irrigated regions occurs during summer when rainfall in the Mediterranean region is very low.Added to this, many of the regions with large-scale irrigation have very dry climates (Lionello, 2012) with the vast proportion of water resources coming from surface water sources.
Increased yield under warmer temperatures and decreased precipitation indicate that in most of the Mediterranean, grain yields are temperature-limited and not water-limited.The spatial distribution of the correlation between climate during the growing season and yield indicates that water is limiting only in very dry regions such as the southern Fertile Crescent, parts of North Africa and coastal regions of the south-eastern Mediterranean (Fig. S7).Increased grain yields under higher temperatures were also found for Mediterranean climate conditions in Western Australia in simulations using the Agricultural Production Systems Simulator (APSIM)-N wheat model to predict the impact of changing temperature, precipitation and CO 2 on yield (van Ittersum et al., 2003;Keating et al., 2003;Ludwig and Asseng, 2006).In the Southern part of the study area (>500mm precipitation), wheat yields were predicted to increase with increasing temperature irrespective of predicted changes in rainfall, whilst in the drier north (<350mm precipitation) rainfall reduction was partially counteracted by increased temperatures (Ludwig and Asseng, 2006).It should be stressed that the response to climate is very heterogeneous throughout the Mediterranean (Fig. S7).
Nonetheless, as we will show, Mediterranean-scale changes are highly relevant at the smaller city-region scale in an integrated network such as the Virtual Water network of the Roman world.

Virtual water redistribution
Rome is the largest importer of VW in our network with imports on average feeding ~460,000 citizens (Fig. 4a).
Egypt is the largest exporter of virtual water, however much of this export is local with large quantities flowing to the densely populated Egyptian cities of Alexandria and Memphis with a proportion also flowing towards Italy (Fig. 4b).The largest flows of VW occur between Eastern and Southern Spain and Rome.There are also large flows between south-eastern Italy and the densely populated region around the Bay of Naples.Other large flows occur along the Turkish Aegean Coast, within the Po Valley and locally in the region around Antioch in present-day southeast Turkey.Although only 28.5% of yield is from irrigated land, VW from irrigated agriculture accounts for 34% of VW flow among nodes.The disproportionately large exports from irrigated land are owing to the location of irrigated cropland close to the coast or along rivers where transport is less costly compared with transport over land.Indeed, all large VW flows originate in areas close to the coast or a large river.Rome has by far the biggest VW demand followed by other large coastal cities such as Alexandria, Ephesus and Antioch (Fig. 4a).
The node degree distribution of the VW redistribution network is shown in Fig. 5a.As with many real-world networks the node degree distribution of our network exhibits a skewed distribution meaning that most nodes are connected to a few edges (low degree nodes) whilst there are a limited number of nodes that are highly connected (high degree or hub nodes) (Konar et al., 2011;Lewis, 2011;Suweis et al., 2011).The correspondence of the node degree distribution to a real-world network gives us confidence that Orbis faithfully captures the network structure of the principal roman trade routes (Scheidel, 2013).Fig. 5b shows the cost to import VW as a function of node degree.
Our analysis indicates that low degree nodes incur the highest import costs in our network (Fig. 5b), consistent with finding that poor infrastructure increases import costs (Limão and Venables, 2001).However in Orbis, lower degree nodes are generally located inland where import costs are also higher owing to the difficulties in transporting large quantities of grain overland by horse and cart compared with ship (Braudel, 1995;Limão and Venables, 2001;Meeks, 2013;Scheidel, 2013).To isolate the effect of node degree from edge cost we simulated VW redistribution with the same network structure but reassigned edge costs and VW values at nodes randomly in each simulation year (Fig. S8a).This analysis demonstrates that import cost is closely related to node degree, independent of the transport costs of edges connected that node.
In a network where costs covary with distance, higher import costs for low degree nodes arise because a node with few transport links has a higher chance of depleting neighbouring nodes compared with a high degree node, assuming equal demand.Once neighbouring nodes are depleted, a VW-poor node must import from further away, thus increasing cost.However, as node degree increases it is less likely that all neighbouring nodes become depleted, which on average will reduce import distance and costs.It is notable that for the highest degree nodes, import costs are higher on average (Fig. 5b).In network theory, highly connected nodes are known as hubs.Hub nodes are mostly located along the Mediterranean coast in our network (Scheidel, 2013).Konar et al. (2011) and Suweis et al. (2011) demonstrated that these hub nodes play a critical role in providing access for poorly connected nodes to the larger VW network.In Orbis, hub nodes are usually ports (for example the port node at Ostia near Rome) or urban centres.
Thus the demand of hub nodes is in reality the sum of demand from many inland nodes or large local populations.
Owing to the high demand levels of these hub nodes they often deplete all their neighbouring VW-rich nodes and must import from further away, thus increasing import costs.
Changes in import costs indicate how stressed our VW network of the Roman World is.For example, if total network cost is 0, then all regions have sufficient local water resources to meet the demands of the local population.
If total network cost > 0 then local water resources in at least one city region are insufficient to meet the local population demand, meaning that VW import is required.To investigate the impact of increased stress on our network, we simulated VW redistribution across a stress gradient based on increases or decreases in population at each VW poor node.We chose to only change populations at VW poor nodes as these are generally representative of urban regions and therefore reflect urban population growth during the late Republican and early Imperial era (Scheidel, 2001).Our analysis indicates that lower degree nodes exhibit a negligible increase in cost as a result of increased demand (Fig. 5c).However, high degree hub nodes exhibit an incremental increase in cost for increasing demand.
In all cases, as demand increases, a VW-poor node must import from further away in the network.For low degree nodes, most of which are inland, the largest costs are involved in bridging the gap to coastal hub nodes.Once a hub node is reached import costs increase relatively slowly owing to the increased number of coastal import routes that can be selected.For high degree nodes, the increased number of import routes that can be selected means that costs 390 begin very low when demand is low and increase incrementally as demand increases and nearby nodes are depleted (Fig. 5c).The outcome of this is that although import costs in poorly connected, inland regions of the network are high, they do not increase substantially for increases in demand.However, for hub nodes that are adapted to low costs, increases in demand can cause substantial increases in import cost.This pattern is only applicable in a network such as our VW network of the Roman World, where lower degree nodes tend to be located inland (Fig. S8b) (Scheidel, 2013), which is also typical of the present-day global trade network for bulk goods (Limão and Venables, 2001).
We find that the total import cost of our water redistribution network is closely linked to climate, in particular temperature (Fig. 6).During warm years, increased yields mean that many regions have sufficient local VW to meet 400 demand so imports are unnecessary.However, even in the case where import is required, total demand will drop with the result that a VW poor node competes with fewer VW poor nodes for increased VW resources.Consequently, nearby surplus nodes are less likely to become depleted and imports occur over shorter average distances.As stated, reconstructions of climate during the Roman period indicate that temperatures were anomalously warm (Chen et al., 2011;Davis et al;., 2003;Ljungqvist, 2010;Wang et al., 2012), creating optimal conditions for the growth of grain.Therefore, the average transport distance of VW in the Empire was likely reduced (Fig. 6) compared with the subsequent, cooler dark ages cold period beginning around 400 AD (Bond et al., 2001;Desprat et al., 2003;McDermott et al., 2001)

Roman Water Resource Management
Taking Rome as an example, our simulations indicate that the majority of its VW was imported from Spain with Sardinia, Southern France and Egypt also contributing substantial quantities (Fig. 4).However, historical sources indicate that Egypt, North Africa and Sicily were the dominant export regions of VW to Rome (Bransbourg, 2012;Erdkamp, 2005).As previously stated, grain yields are underestimated in HYDE for North Africa and Sicily thus Spain supplants these regions as the primary exporters of VW to Rome in the Western Mediterranean in our network.
Additionally, our network solves VW transport along the most efficient routes with VW poor nodes having perfect knowledge of the VW status of the closest VW rich node.Thus import routes are constantly adapted to keep cost to a minimum.However, for the Roman period this is an unrealistic scenario as the efficiency of knowledge transfer varied based on distance, frequency of trade relations etc. (Kessler and Temin, 2007).Probably Romans had little knowledge of yields in remote regions of the Empire until grain ships arrived in port (Rickman, 1980).If the cargo on incoming ships was below expectations then cities risked famine and potential violent unrest (Erdkamp, 2005).
Therefore stable yields in exporting regions were particularly important to the Romans.
Examining the year to year variability in yield we can see that much of the Eastern Empire likely had highly stable yields, in particular Egypt.In the Western Empire North Africa, Sicily and the Po valley exhibited the most stable grain production (Fig. 7).The stability of yields in irrigated regions such as Egypt and the Po Valley was borne out of a year round supply of surface water so that multi-cropping could take advantage of the seasons when temperatures for growth were optimal.Yields from rainfed agriculture were probably most stable in south-western Turkey, the Western Fertile Crescent, North Africa and Sicily.In these regions winter climate is relatively warm compared with Spain, Italy and France and the Adriatic coast (Lionello, 2012).In addition, winter climate was also quite stable owing to the reduced influence of Atlantic Storm tracks compared with the north-western Mediterranean (Lionello, 2012;Xoplaki et al., 2004).Although Spain and France could export large quantities of VW many years, the reliability of yields were much less compared with the aforementioned regions, a disadvantage that was unacceptable in an era of inefficient information transfer (Kessler and Temin, 2007).The high productivity of Spain but low stability in yields is probably why its main exports during the Roman Period were non-staple foods such as olive oil (Blázquez, 1992;Woolf, 1992).
Our analysis highlights that the heterogeneity of the Mediterranean environment was important for providing the Romans with resilience to interannual climate variability.As we mentioned, a widespread use of irrigation as well a warm winter temperatures meant that the Eastern Empire had warmer and more stable winter temperatures compared with the West.Topographical variations also played an important role with grain yields limited by temperature at higher elevations whereas they were water limited at lower elevations and in more arid environments (Fig. S7).It was the Romans ability to link these environmentally heterogeneous regions through trade that provided them with a stable food supply despite the variable climate of the Mediterranean region.However, VW redistribution during the Roman Period also facilitated populations in VW poor regions, in particular urban areas, to overshoot their ecohydrological carrying capacities (Erdkamp, 2005;Garnsey, 1988;Rickman, 1980).The population growth and increased urbanisation during the Late Republican and Early Imperial periods (Bowman and Wilson, 2011;Scheidel, 2001) likely pushed the Empire closer to the limits of available fresh water resources and reduced resilience to climatic variability (D'Odorico et al., 2010;Garnsey, 1988).In addition, our simulations using a cost-distance based network show that increased demand arising from urbanisation caused an increase in average import distance and an associated increase in import costs.It is plausible, that lower water resource redundancy and increased import costs may have been a contributing factor to the third century crisis which followed a period of population growth, urbanisation and trade in the 2 nd century AD (Parker, 1992;Scheidel, 2010).

Present-day implications
Our analysis uncovered a number of important features that have general implications for virtual water trade under spatially and temporally variable environmental conditions.For example, provided there are enough trading regions with temporally heterogeneous yields, virtual water trade increases carrying capacity without an increase in water resource use in any of the trading regions (Fig. 8).Virtual water trade is therefore a highly efficient method of providing resilience to interannual climate variability.However, by increasing the carrying capacity of trading regions as well as allowing VW poor regions to overshoot their local ecohydrological carrying capacities (Fig. 8) (Barnaby, 2009;D'Odorico et al., 2010), virtual water trade promotes population growth and urbanisation (Hubacek et al., 2009).Therefore, the short term resilience that VW trade provides is eroded in the long term as population growth and urbanisation pushes trading societies towards a global carrying capacity.The present-day trend of urbanisation (United Nations, 2014) means that the global society is becoming increasingly dependent on trade to ensure food supplies.As population continues to grow there is less space to adapt to yield perturbations that may arise owing to anthropogenic climate change (D'Odorico et al., 2010).The globalised population is therefore in danger of becoming vulnerable to climate perturbations in the same way that an isolated population is.However, unlike an isolated population, the globalised civilisation is also vulnerable to perturbations in the trade network itself (De Benedictis and Tajoli, 2011;Grubesic et al., 2008).
In our study we used a cost-distance network to investigate VW trade.However, the majority of studies analysing VW trade use socioeconomic trade relations to define the network structure (Hanasaki et al., 2010;Konar et al., 2013;Suweis et al., 2013).Socioeconomic-based trade networks are highly changeable over time with projections of future network structure highly uncertain (Carr et al., 2012a), though evolutionary approaches do show promise (Dalin et al. 2012a).In contrast, cost-distance based networks are much more stable through time.As with the Roman Period, present-day transport costs continue to co-vary with distance, particularly for bulk, staple foods such as grain (Hummels, 2007).Just as in the Roman period, present-day land-based transport of bulk goods is considerably more costly compared with sea-based transport (Limão and Venables, 2001).Indeed, it has been found that trade costs of bulk goods have become increasingly distance sensitive in the latter part of the 20 th and early 21 st century with approximately half of world trade occurring between trade partners less than 3000 kilometres apart (Berthelon and Freund, 2008).The reasons for a stronger relation between cost and distance in recent decades are not straightforward (Berthelon and Freund, 2008), but it is likely that future increases in fuel costs will strengthen the 480 trend further (Curtis, 2009).Therefore, the 'struggle against distance' (Braudel, 1995) which was a characteristic of preindustrial trade remains a central constraint for present-day VW redistribution.
As population growth and urbanisation (United Nations, 2012) continues to reduce water resource redundancy (D'Odorico et al., 2010), our analysis demonstrates that an associated increase in import distance will be unevenly distributed throughout the global VW network, with hub nodes experiencing the greatest increases in import cost.
How such costs will actually be distributed throughout the global trade network are complicated by socioeconomic factors such as protectionism (Carr et al., 2012b;Messerlin, 2011).Therefore, research on VW trade should continue to explore the two-way feedbacks between society and their environment (Sivapalan et al., 2012) because socioeconomic forcings are perhaps the primary driving force of VW trade (Hoekstra and Chapagain, 2008).
However, cost-distance based networks provide an additional avenue for understanding the underlying physical and environmental constraints influencing VW trade.In addition, the stability of cost-distance relations contributes to improving projections as well as identifying the most economical VW trade routes, not just in terms of saving water but also in terms of fossil fuel use.

Conclusion
The question of what brought about the fall of the Roman Empire is one that has occupied Roman scholars for centuries (Gibbon, 1776).However, an equally relevant question is what enabled the Roman civilisation to persist for so long in a region of highly variable climate (Lionello, 2012) and associated variability in agricultural yields upon which their economy and survival depended (Erdkamp, 2005;Garnsey, 1988;Rickman, 1980).Our findings show that the majority of the Mediterranean is temperature-limited for the growth of grain.Given that the height of the Roman civilisation coincided with centuries of anomalously warm climate (Bond et al., 2001;Desprat et al., 2003;Ljungqvist, 2010); conditions for the growth of grain were likely optimal.However, higher frequency climate variability has been demonstrated to have catastrophic impacts for other past civilisations where water resources were not spatially integrated to the extent of the Roman Empire ( de Menocal, 2001;Weiss et al., 1993).The Roman's ability to link heterogeneous environments of the Mediterranean using VW trade meant they could offset deficits in one region with surpluses in another (Kessler and Temin, 2007).In an era before the invention of the combustion engine, the Mediterranean Sea played a critical role because it enabled the spatial integration of the Empire through shipping (Braudel, 1995;Jones, 1986;Scheidel, 2013).The linked-heterogeneity of the Roman Empire undoubtedly increased their resilience to climate variability and contributed to the longevity of their civilisation (Gibbon, 1776;Scheidel, 2013).The importance of VW redistribution in the Mediterranean as a buffer to climate variability is illustrated in the writings of Pliny the Younger (98 AD -117 AD) in Erdkamp ( 2005) 'Even the heavens can never prove so kind as to enrich and favour every land alike.But he [the emperor] can so join East and West by convoys that those people who offer and those who need supplies . . .appreciate . . .having one master to serve'.S1).Virtual water surpluses and deficits are calculated for each region based HYDE reconstructions of Roman Population distribution and grain yields calculated in PCR-GLOBWB.We fix the annual grain consumption at 200 kg per person (D) The virtual water surpluses and deficits are abstracted to the Orbis network and virtual water redistribution is simulated.Deficit nodes import virtual water from surplus nodes along the minimum cost path.Imports continue until the demand of deficit nodes is met.See Fig. S2 for a spatially explicit schematic of our virtual water network of the Roman World.Irrigated regions (E and F) also exhibit increased yields with increasing temperature whereas the impact of precipitation is negligible.The reduced yield under higher precipitation is likely related to decreased temperatures under increased precipitation in most of the Mediterranean and thus lower evapotranspiration (Fig. S6).This indicates that the majority of the Mediterranean is temperature-limited for cereals.For the highest degree nodes, the cost to import is higher than nodes with an intermediate number of links as many of the highly connected nodes in our network are also ports or urban centres with high demand.Therefore nearby nodes are often depleted leading the need to import from further away with an associated increase in cost.(C) For nodes with 1 -4 links, import costs remain high irrespective of the level of demand.However, for nodes with 5 -8 links and 9 -12 links, costs increase under increasing demand.100 percent demand, represents the standard model simulations presented elsewhere in the paper.In an isolated society populations must remain below the climate-forced carrying capacity to avoid famine.In societies with trade, the carrying capacity becomes the average carrying capacity of all trading regions.Thus, where there is a sufficient number of trading regions with heterogeneous environments, carrying capacity is increased without an increase in resource use in any of the trading societies.Trade even allows certain regions to attain carrying capacities well above their local ecohydrological constraints, thus facilitating urbanisation.However, under continued population growth, the resilience to climate variability provided by trade is eroded as populations approach a new global carrying capacity.In addition, urbanisation means that regions become vulnerable to perturbations in the trade network as well as perturbations arising from crop failure.Carrying capacities are smoothed to illustrate the dampening effect of food storage.

Figure 1 .
Figure 1.Schematic of our virtual water network of the Roman world.(A) The overlay of the History Database of the Global Environment (HYDE) cropland reconstructions for 200AD, monthly irrigated and rainfed crop areas (MIRCA) and global land cover characterization (GLCC) forms the land cover attributes in our model.Regions are only assigned irrigated where MIRCA irrigated regions overlap with historical cropland regions.(B) Grain yields in irrigated and rainfed agricultural regions are simulated in PC Raster Global Water Balance model (PCR-GLOBWB) based on 52 years of spatially and temporally variable temperature and precipitation forcing.(C) Regions are defined based on a Theissen polygon operation between network nodes (Fig.S1).Virtual water surpluses and deficits are calculated for each region based HYDE reconstructions of Roman Population distribution and grain yields calculated in PCR-GLOBWB.We fix the annual grain consumption at 200 kg per person (D) The virtual water surpluses and deficits are abstracted to the Orbis network and virtual water redistribution is simulated.Deficit nodes import virtual water from surplus nodes along the minimum cost path.Imports continue until the demand of deficit nodes is met.See Fig.S2for a spatially explicit schematic of our virtual water network of the Roman World.

Figure 2 .
Figure 2. Average cereal yield (Ton per 5' cell).(A) Average cereal yield in tonnes per 5'cell, calculated in PCR-GLOBWB and based on 52 years of climate forcing.The yields from rainfed (B) and irrigated (C) agriculture are shown separately.See Fig. S5 for yield in kg ha -1 .Yields are highest in irrigated regions where year-round supply of surface water allows for multi-cropping, which can take advantage of the seasons when temperatures for growth are optimal.

Figure 3 .
Figure 3. Yield plotted against temperature and precipitation.Total yield (A and B) in the Mediterranean increases with increasing temperature and decreases with increasing precipitation.The trend is strongest in regions where agriculture is rainfed (C and D).Irrigated regions (E and F) also exhibit increased yields with increasing temperature whereas the impact of precipitation is negligible.The reduced yield under higher precipitation is likely related to decreased temperatures under increased precipitation in most of the Mediterranean and thus lower evapotranspiration (Fig.S6).This indicates that the majority of the Mediterranean is temperature-limited for cereals.

Figure 4 .
Figure 4. Virtual water imports and exports.The relative amount of VW imported (A) and exported (B) from each node is illustrated by the size of the nodes, whilst the associated numbers show amount of VW imported or exported in terms of per person population demand at a yearly consumption of 200 kg of grain.The edge colour and thickness indicates the relative volume of VW flow between nodes.The largest flows are between Eastern and Southern Spain and Italy, locally within Egypt, from south-eastern Italy to Western Italy and along the Aegean coast of Turkey.Rome is by far the largest importer of VW, followed by Alexandria and Memphis in Egypt, Ephesus on the West coast of Turkey, Antioch in south-eastern Turkey and Corinth in Greece.

Figure 5 .
Figure 5. Cost to import VW in relation to node degree.(A)The node degree distribution of the virtual water redistribution network.(B) Lower degree nodes generally have higher costs to import VW compared with high degree nodes.For the highest degree nodes, the cost to import is higher than nodes with an intermediate number of links as many of the highly connected nodes in our network are also ports or urban centres with high demand.Therefore nearby nodes are often depleted leading the need to import from further away with an associated increase in cost.(C) For nodes with 1 -4 links, import costs remain high irrespective of the level of demand.However, for nodes with 5 -8 links and 9 -12 links, costs increase under increasing demand.100 percent demand, represents the standard model simulations presented elsewhere in the paper.

Figure 6 .
Figure 6.Cost of imports in relation to temperature.There is a negative relation between the cost to import VW and temperature because yields increase on average in the Mediterranean with increasing temperature.Therefore competition for VW resources reduces and as a result, import distances decrease.

Figure 7 .
Figure 7.The stability of yields over time.The map shows in how many years the total annual yield in each cell remains within 10% of the average yield for the same cell calculated over 52 years of climate forcing.In the Nile Valley, yields remain within 10% of the average yield in all years, meaning that yields are exceptionally stable.Regions of Northern Spain and Northern France are relatively unstable with yields dropping below 10% at least 40 out of 52 years.

Figure 8 .
Figure 8. Conceptual figure illustrating the impact of trade on carrying capacity in a variable environment.Carrying capacities are variable over time owing to the impact of interannual climate variability on yields.In an isolated society populations must remain below the climate-forced carrying capacity to avoid famine.In societies with trade, the carrying capacity becomes the average carrying capacity of all trading regions.Thus, where there is a sufficient number of trading regions with heterogeneous environments, carrying capacity is increased without an increase in resource use in any of the trading societies.Trade even allows certain regions to attain carrying capacities well above their local ecohydrological constraints, thus facilitating urbanisation.However, under continued population growth, the resilience to climate variability provided by trade is eroded as populations approach a new global carrying capacity.In addition, urbanisation means that regions become vulnerable to perturbations in the trade network as well as perturbations arising from crop failure.Carrying capacities are smoothed to illustrate the dampening effect of food storage.