Modelling the risk of deterioration at earthen heritage sites in drylands

The variable climatic and environmental conditions associated with dryland regions can cause rapid erosion to both natural and man‐made earthen structures. Whilst there is a long history of research into the evolution of erosional landforms such as yardangs, little research has investigated how dryland processes influence the erosion of built structures. Earthen heritage sites located in arid and semi‐arid environments experience rapid deterioration caused by exposure to environmental drivers such as wind and rain. Understanding how these environmental drivers interact with each other and cause deterioration to earthen material is vital for successful conservation strategies. To address this need, we present the Vegetation and Sediment TrAnsport model for Heritage Deterioration (ViSTA‐HD) that simulates the risk of polishing, pitting and slurry on earthen heritage in a spatially specific manner. A technical description of the model is provided, and sensitivity and validation tests are reported. The model is then used to simulate the risk of deterioration occurring over centennial timescales at a Suoyang Ancient City, located in semi‐arid northwest China. The modelled risk of deterioration is in good agreement with deterioration patterns found at Suoyang, with the risk of polishing predominantly occurring around the wall edges, areas at risk of pitting echoing the dune formation and the risk of slurry occurring in drape‐like patterns down the wall face. Consequently, ViSTA‐HD is a powerful and versatile model that can be used to help inform our understandings of long‐term interactions between dryland processes and deteriorative impact on earthen structures. © 2020 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd


Introduction
Dryland research has a long history of investigating the climatic and environmental processes that cause the development of natural landforms (Bull, 1977;Wasson and Hyde, 1983;Greeley and Iversen, 1985). The development of erosional landforms, such as yardangs, results from complex interactions between environmental conditions and the parent material (Greeley and Iversen, 1985;Gutiérrez-Elorza et al., 2002). For yardangs, erosion is predominately driven by aeolian processes (abrasion and deflation), with non-aeolian processes (gullying, fluvial processes) influencing the landform's profile (Laity, 2009;Dong et al., 2012;Barchyn and Hugenholtz, 2015). Despite man-made earthen structures in dryland regions being exposed to similar climatic and environmental processes, little research has investigated how dryland processes influence the erosion of earthen structures (Rainer, 2008). These understandings are particularly needed for the hundreds of earthen heritage sites located in dryland regions (WHEAP, 2012), as many are experiencing rapid deterioration due to their low durability to environmental processes (Rainer, 2008;Li et al., 2011).
Wind and rain have previously been found to be the major drivers of deterioration at earthen heritage sites (Wang et al., 2010b;Li et al., 2011;Shao et al., 2013;Du et al., 2017). High-magnitude, low-frequency, wind-driven rain events tend to cause the greatest rate of deterioration, while high-frequency, low-magnitude wind events cause significant deterioration over longer timescales (Shao et al., 2013;Richards et al., 2019). In contrast to many natural landforms that are left to erode, heritage conservators generally aim to minimize the rate of erosion to enable the site to be passed on to future generations. Heritage experts have expressed a preference for earthen heritage conservation strategies to operate over centennial timescales (Richards et al., 2018).
Research into earthen heritage deterioration has typically been undertaken over periods of weeks to years, resulting in the complexity of deterioration processes not being captured over longer time periods (Mao et al., 2015;Cui et al., 2019b;Richards et al., 2019). Deterioration studies undertaken in the laboratory or on test walls have been constrained in their experimental sample size and duration (e.g. Taylor, 1990;Mao et al., 2015;Richards et al., 2019), while field research has been hindered by gaps in (or complete lack of) environmental records for a site. This often results in the classic geomorphological issue of equifinalitywhere two different process histories may feasibly have caused the same deterioration feature (Haines-Young and Petch, 1983). Given these challenges, there is an urgent need for a tool that enables us to: i Understand how the predominant drivers of deterioration (wind, rain and sediment) interact with each other around earthen heritage sites, and what impact these drivers have on the risk of deterioration over a centennial timescale. ii Test the effectiveness of potential conservation strategies over centennial timescales. iii Investigate how future changes to climate may impact the risk of deterioration at earthen heritage sites and test the effectiveness of potential conservation strategies under such climate changes.
Modelling allows us to answer these questions due to its ability to resolve complex processes in dryland environments over multiple spatial and temporal scales (Merritt et al., 2003;Young et al., 2006;Lesschen et al., 2009). Models have commonly been used within geomorphology to investigate interactions and feedbacks between environmental processes (Bagnold, 1941;Schumm and Lichty, 1965;Laity, 2009), but modelling has only rarely been applied to investigate earthen materials and deterioration processes (Nowamooz and Chazallon, 2011;Peng et al., 2015;Bui et al., 2016;Luo et al., 2019). Furthermore, the development of an enviro-heritage model would enable the deterioration risk to be assessed without the need for regular site visits, which, particularly for remote sites, can be expensive and time consuming.
However, modelling interactions between earthen heritage and the environment over a landscape scale can often be computationally expensive. One class of model that has been applied successfully at landscape scale, partly owing to its computational efficiency, is the cellular automaton (CA) model (Nield and Baas, 2008a;Mayaud et al., 2017a). CAs are numerical models based on grid cells that interact with one another through local rules, and have been used to model environmental processes from dunefield evolution (Nield and Baas, 2008a, b;Baas and Nield, 2010;Mayaud et al., 2017a) to dryland vegetation change (Dunkerley, 1997;He et al., 2005;Bailey, 2010;Vega and Montaña, 2011). CA models have been widely applied to simulate processes along the horizontal plane (e.g. Dunkerley, 1997;Kerner et al., 2002;Nield and Baas, 2008a), but vertical-plane processes, such as the deterioration of wall faces, have rarely been investigated using CAs. This is partly a result of vertical-plane processes being generally undertaken over scales small enough to apply computationally expensive models, which are able to capture the complex stresses and strains within vertical materials (Ciantia and Castellanza, 2016;Zhang et al., 2016a;Allen, 2019). However, such computationally expensive solutions cannot be applied to earthen sites, where walls regularly span hundreds of metres.
Given the issues around computational expensiveness of existing vertical models, and the proven ability for CAs to resolve environmental processes horizontally, this paper aims to: i present a CA model that combines a horizontal environmental module with a vertical module that determines the risk of deterioration across a wall face; ii verify and validate the model performance using data collected from an earth heritage site.
We use the Vegetation and Sand TrAnsport model (ViSTA), originally designed to investigate sand dune reactivation in the Kalahari (Mayaud et al., 2017a,b), as a basis for developing the new Vegetation and Sand TrAnsport model for Heritage Deterioration (ViSTA-HD). We combine dryland geomorphological concepts with data collected at the earthen heritage site of Suoyang Ancient City (northwest China) to build deterioration risk functions into ViSTA-HD. We present a range of sensitivity and validation exercises to ensure the reliability and accuracy of our model. The ViSTA-HD model we present here is a first attempt at testing the feasibility of coupling environmental processes with deterioration risk of earthen heritage. ViSTA-HD has been designed as a tool to inform conservation strategies over decadal to century timescales, allowing users to evaluate potential consequences of climate change for earthen heritage. Although the present form of ViSTA-HD is based on environmental and deterioration data from a specific earthen heritage site, its functional versatility allows it to be adapted and applied to other earthen sites in dryland environments around the world.

Data Collection
The behaviour of the model is driven by forcing conditions (e.g. initial windspeed) and parameters used to represent various physical processes. As processes are not resolved explicitly, we use parameterizations to account for the effects of the processes, chosen with the aid of empirical data. We first report on a series of in-situ experiments designed to provide such data. We collected data on environmental conditions that can cause deterioration to earthen material (e.g. wind velocity, wind flow round vegetation and sediment transport), as well as on the type and extent of deterioration features present on earthen heritage.

Field site
Field experiments were undertaken at Suoyang Ancient City (锁 阳城), hereafter known as Suoyang, during a 2-week fieldwork period in August 2018. Suoyang is located in the semi-arid area of Gansu Province, China, and dates back to the Han (206 BC-220 AD ) and Tang (618-907 AD ) dynasties (UNESCO, 2014). It is an archaeological site built out of rammed earth that was aban-doned~400 years ago (Xie et al., 2007). The site was enlisted as part of the Silk Roads World Heritage Site in 2014 and exhibits a variety of deterioration features caused by wind, sediment transport and rain (Richards et al., 2019). This has resulted in pitting, polishing and slurries occurring on the wall faces, as well as the formation of larger features on the earthen walls such as cracks, gullies, slumps and collapses ( Figure 1).
Suoyang's remote location and harsh continental climate mean that it has remained understudied until recently. Since 2015, an experimental test area has been constructed at the site, including a set of four large test walls each measuring 3 m long, 2 m high and 1 m wide. As a result, the site of Suoyang now provides a unique opportunity to investigate environmental and deteriorative processes representative of many other remote sites in this climatic region (e.g. Wang et al., 2010a;Li et al., 2011).

Wind
Wind flow was measured around one of the four test walls at Suoyang. A transect of cup anemometers (height 0.5 m) was set up perpendicular to the test walls (similar to the setup of Mayaud et al., 2016a), with anemometers located at distances of 0.5h, 3h, 5h, 10h and 20h (where h is the wall height) upwind and downwind of the wall (Figure 2a). Two additional anemometers where placed on each side of the wall (one 2402 J. RICHARDS ET AL. between the two walls and the other on the open side of one of the walls). Wind velocity was recorded every minute for 7 days. Measurements were filtered according to the approaching wind direction, and average wind velocities were normalized to measurements at 20h upwind (Van Boxel et al., 2004). Two wind towers with anemometers 0.2, 0.5, 1, 1.5 and 2 m above the ground were placed at 1h either side of the test wall ( Figure 2a; Tsoar, 1983). Measurements were taken every minute for 18 days and filtered according to the approaching wind direction (Van Boxel et al., 2004).
To capture the turbulent nature of wind flow around the walls, flares loaded with theatrical smoke were deployed upwind of the test walls when the prevailing wind was perpendicular to the test wall (Piscioneri et al., 2019). The flares lasted approximately 3 min. Smoke flow was video-recorded, and stills from the videos used to capture turbulent flow patterns in front and around the wall face (Livingstone, 1986;Piscioneri et al., 2019). Macro-scale smoke patterns were used to visualize sediment transport pathways around the walls.

Sediment Transport
Two types of sediment traps were deployed for 1 month between June and July 2018 (Gares et al., 1996). Sediment traps, similar in design to a wedge-shaped sediment trap (0.9 m tall with 45 × 0.02 m increment gaps), were deployed along the same transect as the anemometers facing the prevailing wind direction (Figure 2a; Nickling and McKenna Neuman, 1997). In addition, a circular sediment trap (0.5 m tall) with 16 directional openings (0.02 m wide) was deployed to the east of the site (Greenley et al., 1996). The circular trap was located away from the walls, with a clear fetch in all directions. After 1 month, the mass of sediment captured was recorded for each 0.02 m section in the wedge-shaped sediment trap and for each directional opening in the circular trap. All samples were subjected to grain size analysis (Arens et al., 2002) using a Winner 2,308 Dry Wet Particle Size Analyzer.
To gain a broad-scale understanding of sediment movement around the historic city, we used a laser distance measurer to record the dimensions (height, width and length) of small dunes that had formed against the inner walls of Suoyang.
The sediment data were used to inform the sediment transport within the model space and how the volume of sediment varied with wall height.

Vegetation
Vegetation elements in dryland environments impact wind flow patterns in a variety of ways (Wolfe and Nickling, 1993;Mayaud et al., 2016b). The porosity of vegetation is one of the primary factors determining the extent to which vegetation slows down wind in its lee and impacts flow recovery (Plate, 1971;Wolfe and Nickling, 1993;Torita and Satou, 2007;Mayaud et al., 2016a). In this study, optical porosity was determined using a digital photogrammetry method adapted from Mayaud et al. (2017a), where photographs of each vegetation element were converted into greyscale images. Black (plant) and white (open) areas were selected using a visually determined threshold value. The mean optical porosity was equal to the percentage of white pixels.
Initial optical porosity tests showed that the structure of the vegetation at Suoyang was similar to vegetation in other semi-arid areas (e.g. in the Kalahari; Mayaud et al., 2017a), although some vegetation elements were notably denser. Vegetation elements at Suoyang with similar porosities to those studied by Mayaud et al. (2016a) were assumed to have the same impact on wind flow. However, to test the impact of higher-density vegetation on wind velocity, wind velocity was measured around five high-density vegetation elements (with porosities ranging from 2 to 23%) using cup anemometers deployed in the same manner as around the test walls ( Figure 2b; Mayaud et al., 2016a). Wind measurements were taken by each cup anemometer every 30 s for between 40 and 50 min for each vegetation element to ensure a sufficient sampling period was captured after the measurements were filtered according to the approaching wind direction (Van Boxel et al., 2004). Filtered average wind velocities were normalized to measurements at 20h upwind (Richards and Hoxey, 2012;Mayaud et al., 2016a).
Vegetation change at Suoyang was assessed using image analysis (Alves et al., 2003;Youhao et al., 2007;Duguay et al., 2014;Rehman et al., 2019). Satellite images from 2006 and 2010 were compared to drone photos taken of the site in 2018. Vegetation was manually identified and determined as a percentage of the total area by digitally calculating the ratio of white pixels (no vegetation) to black pixels (vegetation) (Kenney, 1987).
Vegetation data were used to inform the wind recovery distances around vegetation elements within the model, and whether vegetation coverage should remain static in the model runs reported below.

Average climatic conditions
Average wind velocity was recorded every 10 min by an automatic weather station at Suoyang, intermittently between 2015 and 2018 (Richards et al., 2019). Wind measurements taken between September 2015 and August 2016 were used to construct probability density functions which, for each model iteration, provide wind velocity values for wind entering the model domain, and then as a reference throughout the domain during calculation of wind flow interruption ( Figure S6).
Rainfall data collected by the automatic weather station between March and December 2018 were used in conjunction with historical climate data for the area provided by the Dunhuang Academy (unpublished) to inform the rainfall parameters in ViSTA-HD ( Figure S7).

Deterioration data
A semi-quantitative deterioration assessment was undertaken at Suoyang, which included both the type and extent of deterioration present on the earthen walls (Zhao et al., 1994;Turkington et al., 2003). Fifteen survey sites were selected around the inner-city wall to capture deterioration across the site and in a range of microclimates (e.g. on both highly exposed and sheltered walls). At each survey site the presence and extent of deteriorationincluding polishing, pitting and slurrywas noted (Richards et al., 2019), and photographs were taken for further desk-based analysis (Viles and Goudie, 2007;Shao et al., 2013). Richards et al. (under review) provide more information on the deterioration assessment. Outputs from ViSTA-HD under a range of climatic conditions were compared to the deterioration assessment to test ViSTA-HD's ability to capture the risk of earthen deterioration under multiple different climate scenarios.

Model Description
ViSTA-HD comprises two modules: (i) an environmental module that spatially resolves environmental processes occurring around earthen heritage; and (ii) a deterioration module that resolves the risk of deterioration occurring across a wall face ( Figure 3). For each iteration, ViSTA-HD determines the environmental conditions at the upwind face of the wall. These conditions are inputted to the deterioration module to calculate the deterioration risk. ViSTA-HD is designed to simulate environmental conditions occurring along the upwind wall face, because this zone displays the most severe and complex deterioration patterns at Suoyang (Richards et al., 2019). For the purposes of this study, we focus on the three most prevalent forms of deterioration at Suoyangpolishing, pitting and slurry (Richards et al., 2019)but our versatile modelling framework allows other environmentally driven deterioration processes to be added if required.
ViSTA-HD is implemented in the Python® programming language (Python 3.7.0 64 bit). It uses code written by the authors that is openly available from the Oxford University Research Archive (DOI: 10.5287/bodleian:nr86P0OjY). 2404 J. RICHARDS ET AL.

Environmental module
The environmental module is based on the ViSTA model (Mayaud et al., 2017a), which was designed to explore vegetation and dune dynamics in the semi-arid Kalahari Desert. The ViSTA-HD environmental module builds on ViSTA by adding parameterizations of environmental processes occurring at Suoyang, including solid earthen walls in the model space.
ViSTA-HD's environmental module explicitly resolves wind flow and sediment transport interactions with elements such as vegetation and earthen walls across a horizontal model domain comprised of uniform grid cells. In this paper, we outline the relevant functions that have been adapted for the new ViSTA-HD model; see Mayaud et al. (2017a) for a detailed description of the original ViSTA model.

Vegetation
The wind flow in the lee of elements was parameterized in ViSTA as where u surf is the surface wind velocity (m s À1 ) at a given distance downwind (x m) of u ref , an undisturbed reference flow. u 0 is the mean normalized minimum wind velocity at 1h in the lee of the element, x h is the normalized distance downwind from the element and b is the fitted recovery distance (Mayaud et al., 2017a). As u 0 and b were derived from experimental data and are dependent on porosity, we parameterize them in this study using a combination of datasets from Suoyang and Mayaud et al. (2017a) (see Figure S2). The mean normalized minimum wind velocity (u 0 ) was parameterized as where p is the porosity of the vegetation element. The rate of wind velocity recovery behind a vegetation element (b) is parameterized using a threshold relationship, where b¼ 0:32 Field studies have found small increases in wind speed around the edges of vegetation elements (Mayaud et al., 2016b). However, given the limited magnitude of the acceleration and the small width of the acceleration zone around the plant compared to the model resolution, acceleration edge effects around vegetation elements have not been included in this version. See Section S1.1.2 in the online Supporting Information for more details.

Solid walls
Wind velocity measurements taken around the test walls showed that wind velocity decelerates from roughly 3h upwind of the wall (where h is the height of the wall) owing to the formation of a turbulent zone in front of the wall (Tsoar, 1983). Wind velocity accelerated around the edge of the walls, and the magnitude of acceleration was greater in the gap between two test walls (Baskaran and Kashef, 1996;Blocken and Carmeliet, 2004b).
Based on our experimental data, the acceleration up to 1 m around the edge of the wall (W edge ) is parameterized in the model as where W À1h is the windspeed at 1h upwind of the test wall and where θ is the wind direction in degrees perpendicular to the wall face ( Figure S3).In the lee of the test walls, wind recovery is determined using the same parameterization as applied to vegetation (i.e. the wall is considered to be an element with 0% porosity). Research on valleys has shown that incidence angle can impact wind recovery distances (Wiggs et al., 2002;Garvey et al., 2005), but given the small scale of walls, we keep the recovery distance constant following Mayaud et al. (2017a).
The flare experiments and scour markings in front of the walls showed that wind flow interaction with a solid object caused sediment transport pathways to be deflected along the wall face ( Figure S4). Therefore, sediment transport pathways in the model are linear unless they intercept an earthen wall, in which case the pathway gets deflected along the wall face (see Figure 4).
Temporal scale ViSTA-HD runs for a number of iterations set by the user. The period of time represented by each iteration, and thus the 'real time' simulated in a run, is dependent on the sediment flux rate, which is also set by the user. To calculate the duration of each iteration, we undertook a verification test comparing the total sediment flux in each modelled iteration to the sediment flux at Suoyang estimated using two techniques: (i) calculated using Dong et al.'s (2003) sediment transport model, which quantifies the amount of sediment flux expected from an aeolian system under ideal conditions; and (ii) calculated from the mass of sediment collected in the field using a circular sediment trap.
The sediment flux in ViSTA-HD was determined by running the environmental module with environmental and climatic conditions representative of those at Suoyang (Figures S6 and  S7). ViSTA-HD was run 10 times for 200 iterations. The median sediment flux per iteration was 589 kg m À1 . Dong et al.'s (2003) model, reformulated by Bailey and Thomas (2014), was used to calculate the total sediment flux using wind speed and sediment particle size data collected at Suoyang as inputs. Sediment flux was determined by where q(u) is the sediment flux with a wind speed of u (m s À1 ), c 1 and c 2 are fitted constants, d is the grain size (95 μm), D is the reference grain size (250 μm), g is the acceleration due to gravity, p is the air density and v T is the velocity for entrainment  Figure S1), each iteration was set to represent 3 months.

Process interactions
Feedbacks between rainfall events and sediment entrainment were not included in ViSTA-HD for three reasons. First, Suoyang receives very low precipitation, with an averagẽ 1 mm per month in winter and~9 mm per month in summer. Second, rainfall events occur infrequently: rain fell on a single day in SON 2018, and on 12 days in JJA 2018. Third, the  surface sediment layer (top 0-10 cm) usually dries within 24 h after rain events in semi-arid environments (Ravi et al., 2006;Okin, 2016). In fact, during our fieldwork, 10 mm of rainfall fell in a 36-h storm but the surface sediment had totally dried by the following day. Therefore, the period of inactive sediment transport is considered to be negligible compared to the seasonal scale (3 months) at which the environmental module is run.
Image analysis of vegetation cover at Suoyang revealed very little change in vegetation cover between 2006 and 2018 (range of 52-54% cover). Consequently, rainfall-vegetation feedbacks were also not included in our experiments, and vegetation cover remained constant throughout each model run. Running ViSTA-HD at other earthen heritage sites may require the inclusion of vegetation growth, which is easily 'turned on' in the environmental module.

Deterioration module
The authors added the new deterioration module to the modified ViSTA model as a vertical CA that simulates the deterioration risk of polishing, pitting and slurrythe three major deterioration features seen at Suoyang (Richards et al., 2019) and other earthen heritage sites in northwest China, such as Jiaohe Ancient City, Xinjiang Province (Shao et al., 2013) and Yumen Fort, Gansu Province (Wang et al., 2010a). Polishing is caused by the removal of loose grains of material, resulting in a smooth surface. Pitting occurs when environmental processes loosen and deflate grains from the parent material, resulting in concave pits forming on the surface (Mao et al., 2015) and slurry is caused by the deposition of runoff-borne sediment down the wall face, resulting in the formation of a harder crust (Cui et al., 2019a).
In each iteration, the probability of polishing, pitting and slurry occurring in a cell is calculated according to 1 The impact of the environmental drivers in that given cell (environmental risk; Figure 5a). 2 The risk of deterioration in neighbouring cells (spatial risk; Figure 5b). 3 The risk of deterioration that previously occurred in the given cell (temporal risk; Figure 5c).
The calculated probability is then used to determine the deterioration risk in the given cell (see later). The cells in the deterioration module are set at a finer spatial resolution than those in the environmental module (0.25 and 1 m 2 , respectively), reflecting the difference in scale between processes operating over a landscape compared to processes causing a risk of deterioration across a wall face (Shao et al., 2013;Cui et al., 2019a;Richards et al., 2019).

Environmental risk
The environmental risk (E risk ) is calculated using risk matrices as shown in Figure 5a. The environmental risk of polishing and pitting is determined by wind velocity and sediment transport, whereas the environmental risk of slurry is determined by rainfall and wind velocity.
Wind. Wind velocity is resolved at 50 cm above the surface in the environmental module. Therefore, the wind velocity is inputted to the deterioration module at 50 cm above the surface ( Figure S5).
The increase in wind velocity up the wall face was parameterized using the data from the anemometry towers located where W h is the wind velocity (m s À1 ) at a given height up the wall, W 50 is the inputted wind velocity from the environmental module at 50 cm above the surface and h is the height above the surface (cm) ( Figure S6a).
Sediment. All sediment deposited in front of a wall face (as determined by the environmental module) is assumed to have interacted with the wall. Our sediment trap data suggests that the volume of transported sediment decreases with height ( Figure S6b), with 67% occurring in the first 40 cm above the surface. This compares well with previous studies, which showed that saltation occurs up to 40 cm above the surface (Dupont et al., 2013). Thus, we set two-thirds of the sediment to interact with the wall, within the height limit of the saltation zone. (8) where S h is the volume of sediment hitting the wall surface at a given height (cm), S total is the total volume of sediment deposited in front of the wall per iteration, S 0 À 40 is the total volume of sediment that interacted with the wall in the first 40 cm and h exp is the height of the exposed wall above any local accumulation of sand.
Rainfall. The impact of rainfall on earthen heritage is dependent on rainfall volume and wind velocity. Below 2 m s À1 , rain is considered as only hitting the wall top in accordance with BS EN ISO 15927-3:2009. At, and above, 2 m s À1 , the wind provides enough horizontal momentum for wind-driven rain to occur. Whilst research on buildings has shown that the base of walls are more protected than the tops (Blocken and Carmeliet, 2004a), many earthen heritage structures are shorter than western-style buildings and do not have typical building attributes, such as roofs, which affect the intensity of wind-driven rain across a wall. Furthermore, surveys of earthen heritage after storm events have shown moisture impacts across the entire wall face (see Section 2.2 in the online Supporting Information). Therefore, wind-driven rain is modelled as evenly affecting the wall face.

Spatial risk
The spatial risk of deterioration captures spatial vulnerability to deterioration. As the material properties of the earthen wall (e.g. variations in clay content, grain size, etc.) are not explicitly determined in this model, the spatial risk captures the principle that if neighbouring cells are deteriorated, a given cell is more likely to be vulnerable to deterioration. It is modelled as a linear function: where S risk is the spatial risk, det is the type of deterioration, X is 2408 J. RICHARDS ET AL.
a scaling factor, Nc det is the number of deteriorated neighbouring cells and Nc total is the total number of neighbouring cells. For polishing and pitting, neighbouring cells are considered to be one cell on each side of the polled cell (including diagonals; Figure 5b). For slurry, neighbouring cells are considered to be one cell above and to the side of the polled cell (including upper diagonals; Figure 5b) due to the role of gravity in bringing moisture down the wall face (Zhang et al., 2016b).

Temporal risk
The temporal risk (T risk ) is determined by the deterioration that occurred in the polled cell on the previous iteration. Previous deterioration processes can prime wall surfaces, reducing and enhancing their susceptibility to future deterioration (Figure 5c; Cui et al., 2019a). The temporal risk values were determined by results from previous studies (e.g. Du et al., 2017;Cui et al., 2019a) and on-site visual deterioration assessments at Suoyang. They are modelled as a linear function, as shown in Figure 5c.

Deterioration risk
The risk of pitting, polishing and slurry is calculated separately in each cell by multiplying the environmental, spatial and temporal risks to create three different risk probabilities, P(det) ( Figure 5d): The cell is then assigned a value of either 0 or 1, based on r¼random number between 0 and 1 cell value¼ In each cell, for each deterioration type, these binary outcomes are collated and averaged over all model iterations (n): The authors would like to strongly emphasize that a value of '1' does not mean that deterioration will definitely occur and '0' does not mean that deterioration will definitely not occur. These values are merely designed to give an indication of relative risk. For the purposes of this study, we consider the risk of deterioration as the potential for deterioration to occur in a specific special areanot the severity of material loss.

Sensitivity Testing of ViSTA-HD
We present the results from the sensitivity testing undertaken on ViSTA-HD. Additional model verification tests to show that the assumptions and parameterizations used when building ViSTA-HD are suitable for the purpose of investigating earthen heritage deterioration are presented in Section S2 of the online Supporting Information. Further verification tests of the environmental module can be found in the Supporting Information of Mayaud et al. (2017a). In all sensitivity tests, to balance computational expense with capturing smaller-scale deterioration features, a resolution of 0.25 m 2 was chosen for the deterioration module (see Section S2.2.3 in the online Supporting Information).
The sensitivity of the deterioration module was tested by varying the magnitude of the environmental, spatial and temporal risks ( Figure 6). All combinations of environmental, spatial and temporal risk scenarios were cross-tested, resulting in 27 outputs (Figure 7). Each scenario was run with the same set of starting conditions, representative of those found at Suoyang, repeated 50 times until the mean and standard deviation had equilibrated (see Sections S2.1 and S2.2.2 in the online Supporting Information), and then averaged. Zones of the wall face with high levels of agreement across model runs (i.e. standard deviation < 5% of the mean risk of deterioration) are hatched, with the variability in deterioration risk generated in Equation 11. Areas which are exposed to very high or very low risks of deterioration have the least variability in deterioration risk, whilst areas at a medium risk of deterioration have greater variability.
The sensitivity tests resulted in significantly different outcomes when the environmental conditions were run under the uniform scenario ( Figure 7a) compared to the low ( Figure 7b) and high scenarios (Figure 7c). The uniform scenario resulted in high risk probabilities with equally distributed risk across the wall face, which is unrepresentative of deterioration patterns observed at Suoyang. Both the low and high environmental risk scenarios showed similar patterns of deterioration across the wall face, with the low scenario having a notably lower distribution of deterioration risk (Figure 8a). The high deterioration scenario showed greater consensus between model runs and captured a wider range of risks, which were more representative of the deterioration features observed at Suoyang.
There were also notable differences between the results from the spatial and temporal risk scenarios. The 'no impact' spatial and temporal risk scenarios resulted in the lowest deterioration risk, with the high scenarios having the greatest increase in deterioration risk (Figures 7a-c). As both spatial and temporal risks are multiplication factors applied to the environmental risk, the impact of the spatial and temporal risk increases with the initial environmental risk. Therefore, if the spatial and temporal scales are set too high, the scale may 2409 MODELLING RISK OF EARTHEN HERITAGE DETERIORATION become saturated in highly deteriorative environments. This is evident from the peak in the distribution around 1 for high spatial risk (Figure 8b). For high temporal scale scenarios, the distribution of mean deterioration risks is more widely spread between 0.75 and 1, indicating that saturation has not occurred (Figure 8c).
Given the results from our sensitivity analyses, the standard conditions for the deterioration module are set to those corresponding to a high environmental risk, a low spatial risk and a high temporal risk ( Figure 5).
Future efforts in the field should be directed at further quantifying the deterioration risk (i.e. more precisely defining what FIGURE 6. The probabilities used in the sensitivity testing for (a) environmental risk scenarios uniform, low and high; (b) spatial risk scenarios none (no impact), low and high; and (c) temporal risk scenarios none (no impact), low and high. [Colour figure can be viewed at wileyonlinelibrary.com] FIGURE 7. The mean risk of deterioration after 50 model runs when inputting (a) uniform, (b) low and (c) high environmental conditions subjected to varying spatial and temporal risks (ER, SR and TR, respectivelysee Figure 6 for definitions of uniform, none, low and high). Hatching indicates areas of the wall where the standard deviation is <5% of the mean. Note the different scales used for ER uniform, compared to ER low and high.
environmental, spatial and temporal conditions lead to specific deterioration features). Such data would enable a quantification of the impact of uncertainty in these parameters on the final model outcomes. In the meantime, we chose parameters for the environmental, spatial and temporal risk that yield sensible outcomes in relation to field conditions. We note that our sensitivity analysis shows that model outcome is not sensitive to these precise values, such that any uncertainty in our final results will not be strongly impacted by uncertainty in these parameters.

Model Validation
Deterioration patterns at Suoyang are used to validate: (i) the risk of individual deterioration features; and (ii) the overall risk of deterioration under five different environmental conditions occurring on earthen heritage.

Formation of deterioration features
When inputted with environmental conditions representative of those at Suoyang, ViSTA-HD produced deterioration patterns of polishing, pitting and slurry observed at our field site ( Figure 9). Polishing dominates at wall tops and edges (Figure 9b-i), pitting follows the development of dunes up the wall face (Figure 9c-i) and slurry drapes down the wall face (Figure 9d-i).
In this simulation, given the higher modelled risk of slurry compared to polishing and pitting, slurry patterns dominate in the overall deterioration, except for at the wall edges where polishing plays an important role (Figure 9e).

Simulating deterioration in different environments
ViSTA-HD is able to capture the deterioration risk on walls exposed to different microclimates around Suoyang ( Figure 10). Five scenarios were tested to simulate areas: i around the on-site weather station exposed to seasonal variations in wind velocity and rainfall; ii sheltered from the prevailing wind and rain; iii exposed to the prevailing wind and sheltered from the rain; iv sheltered from the prevailing wind and exposed to the rain; and v exposed to the prevailing wind and rain.
See Figure 10 for definitions of sheltered/exposed wind/rain.
Under conditions similar to weather station records, deterioration was generally evenly distributed across the wall face, with some light drape features forming (Figure 10a). This is observed at Suoyang, where many of the walls exhibit some minor deterioration across the wall face, with notable features occurring in specific areas such as around weaknesses or cracks in the material (Figure 10a-i). Under sheltered wind and rain conditions, there was minimal risk of deterioration except for exposure of the top corners to polishing (Figure 10b). With increased wind velocities, the areas at risk of polishing increase but remain predominantly located near the wall tops. This is highly reflective of the deterioration observed on walls sheltered from wind and rain at Suoyang that exhibit small amounts of pitting and polishing along their face (Figure 10b-i). Under exposed wind and sheltered rain conditions, the risk of polishing increases notably on the upper half and on the corners of the wall (Figure 10c). Whilst most areas at Suoyang exposed to high winds are also exposed to wind-driven rain, polishing and pitting zones occur at the top of walls that are more sheltered from the rain (Figure 10c-i).
When subjected to sheltered wind and exposed rainfall conditions, drapes of slurry emerge along the wall face (Figure 10d). This pattern corresponds well to thick slurry flows observed in sheltered areas of Suoyang's walls (Figure 10d-i). In exposed wind and rain conditions, the impact of wind-driven rain is seen with greater overall risks of deterioration ( Figure 10e). However, in this case, the drapes of slurry are less distinct than in Figure 10d; instead, a more homogeneous sheet of slurry forms across the wall face. This particular feature is seen at Suoyang along walls that are exposed to high wind and rain (Figure 10e-i).

Applications of ViSTA-HD
ViSTA-HD is a robust model that can simulate environmentally driven deterioration risk across earthen walls. It can simulate a range of deterioration risk patterns representative of the deterioration patterns at Suoyang. ViSTA-HD enables environmental drivers of deterioration to be assessed over a range of temporal and spatial scales, which is not always feasible in laboratory or field tests. The model also allows practitioners to test the outcomes of different nature-based conservation strategies over decadal and centennial timescales.
ViSTA-HD has the potential to be adapted and applied to other earthen sites in dryland environments, with the addition of site-specific data such as vegetation coverage and sediment flux rates, and the inclusion of other environmentally driven deterioration features if prevalent. ViSTA-HD can thus act as a decision-making tool to help conservators and site managers assess trade-offs between possible interventions, enabling a holistic assessment of strategies that are likely to be most effective today and into the future.
ViSTA-HD also acts as a tool for studying the evolution of natural landforms made out of earthen or semi-lithified materials in dryland environments. For example, ViSTA-HD could be adapted to explore the formation of erosional features such as yardangs under present as well as past and future environmental conditions (Ward and Greeley, 1984;Dong et al., 2012).

Conclusion
Understanding the complex interactions between deterioration processes at earthen heritage sites requires tools that explicitly account for temporal and spatial scaling. With this in mind, we developed a new enviro-heritage model, ViSTA-HD, which can spatially resolve the risk of deterioration resulting from wind, sediment transport and rain across an earthen wall. In its current form, the model produces realistic examples of three forms of deteriorationpolishing, pitting and slurrythat arise from exposure to different environmental conditions over space and time. This could be expanded to include other forms of deterioration, if required. Consequently, ViSTA-HD is a powerful and versatile model that can inform understandings of long-term earthen site deterioration and address the researchconservation scale gap currently facing many earthen heritage sites. By reducing the burden on researchers to carry out fieldwork in remote and difficult conditions, models such as ViSTA-HD are invaluable tools for testing the potential outcomes of conservation strategies, and for predicting the likely impacts of future climate and land use change on earthen heritage.
The ViSTA-HD model is now ripe for further testing and validation. In future work, we envisage ViSTA-HD being used to investigate deterioration risk at other earthen heritage sites in semi-arid and arid environments, and having the potential to be adapted for understanding the development of natural erosional landforms.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. Wind velocity recorded at Suoyang between September 2015 and August 2016, grouped into four seasons. Figure S2. The relationship between porosity and (a) the mean normalized minimum wind velocity and (b) the b coefficient using the joint Kalahari-Suoyang dataset. Figure S3. (a) The relationship between wind direction and relative wind velocity at two locations shown in (b). The relative wind velocity is calculated according to the wind velocity at 1h upwind of the wall (where h is the height of the wall). Figure S4. A photograph of the flare experiment showing the deflection of the wind flow around the test wall. Figure S5. (a) The interaction of the environmental and deterioration modules at the wall face: x 1 and x 2 represent environmental conditions resolved in the environmental module. The solid black arrows show interpolation between inputs. Dashed vertical arrows indicate the extrapolation of conditions above and below the environmental input and the dashed horizontal arrows indicate the interpolation between calculated conditions across the wall face. (b) An example of wind velocity inputs from the environmental module into the deterioration module. Figure S6. Change in (a) wind velocity (relative to the wind velocity 0.5 m above the ground) and (b) total volume of transported sediment, with increasing height from the ground. The dashed black line in (b) indicates the determined height for sediment transport occurring within and above the saltation zone. Figure S7. Wind roses comparing the wind regime above 4 m s À1 at Suoyang to the wind velocities inputted into ViSTA-HD. Figure S8. (a) Total rainfall, (b) quantity of rainfall in each storm event and (c) the number of rain events per iteration, over an annual time period and for each season (DJF, MAM, JJA, SON) individually. The dashed black lines in the first panel show the mean annual rainfall and the solid black lines shown the minimum and maximum annual rainfall. Figure S9. A comparison of the modelled wind flow around three walls of different heights (3, 5 and 8 m) with data collected from Suoyang: (a) location and heights of the three modelled walls; (b) wind velocity across the model space caused by a wind blowing from left to right; (c) transects of wind velocity to the side of the wall where acceleration occurs and through the centre of the wall; (d) comparison of modelled wind velocities (grey) to wind velocities collected at Suoyang (black). The distance from the wall is reported in h (where h is the height of the wall), to enable comparison between model and field data. Error bars show one standard deviation from the mean. Figure S10. The development of dunes in front of two walls (measuring 50 and 10 m in length and 6 m in height), over 300 iterations. Figure S11. Comparison of modelled and real dune morphologies on the upwind side of walls. The range of dune morphologies are shown using histograms for the modelled dunes and triangles for the measurements taken in the field. Four model configurations were used to compare (a) the height of exposed walls and (c) the dune length. An example of each model configuration after 500 iterations is shown in (b). The model setup is shown in (d). Figure S12. Comparison of dune morphologies between models starting with a uniform (flat) sand height and a pre-conditioned surface for (a) dune length and (b) exposed wall, where h is the height of the wall. Figure S13. The impact of rain on (a) test wall 24 h after a storm event and (b) historic walls. Figure S14. The running mean for 50 repeat runs of the (a) mean and standard deviation of the risk of deterioration across two wall faces and (b) standard deviation of the risk of deterioration in randomly selected individual cells. The approximate locations of the individual cells are shown on wall 1. Figure S15. Overall risk of deterioration modelled at 2 cm resolution on walls representing walls in locations (a) sheltered from the prevailing wind and sheltered from rainfall, (b) exposed to the prevailing wind and sheltered from rainfall, (c) sheltered from the prevailing wind and exposed to rainfall and (d) exposed to the prevailing wind and exposed to rainfall. Here, sheltered locations are modelled with a mean wind velocity of 5 m s À1 and the median rainfall event being 1 mm with an interquartile range of 1.75 mm. Exposed locations modelled with an average wind velocity of 12 m s À1 and the median rainfall event being 3 mm with an interquartile range of 8 mm.