Application of an entropy maximizing and dynamics model for understanding settlement structure: the Khabur Triangle in the Middle Bronze and Iron Ages q

We present a spatial interaction entropy maximizing and structural dynamics model of settlements from the Middle Bronze Age (MBA) and Iron Ages (IA) in the Khabur Triangle (KT) region within Syria. The model addresses factors that make locations attractive for trade and settlement, affecting settlement growth and change. We explore why some sites become relatively major settlements, while others diminish in the periods discussed. We assess how political and geographic constraints affect regional settlement transformations, while also accounting for uncertainty in the archaeological data. Model outputs indicate how the MBA settlement pattern contrasts from the IA for the same region when different factors affecting settlement size importance, facility of movement, and exogenous site interactions are studied. The results suggest the importance of political and historical factors in these periods and also demonstrate the value of a quantitative model in explaining emergent settlement size distributions across landscapes affected by different socio-environmental causal elements.


Introduction
In the Middle Bronze Age (2000e1600 BC; MBA), large settlements arose in northeast Syria in the region of the Khabur Triangle (KT). The distribution of settlement sizes in this region was relatively broad; numerous small and medium sized sites arose, while few sites became very large. These structures reflect the rise of kingdoms that briefly integrated this region within larger states, while at other times the region was politically fragmented. Roughly one millennium later in the KT, the political and settlement picture in the Iron Age (1200e600 BC; IA) were very different when the large territorial Neo-Assyrian state dominated the landscape. Although during this long period there are likely key settlement size variations, very few settlements reached sizes evident in the MBA and many sites were small.
It is not well understood how such settlement size structures developed based on inter-and intra-regional interactions and socio-environmental factors. While there are historical sources that explain the political situation for the region, we are missing vital information. Based on this, a methodology is needed that can explain such urban distribution structures, notwithstanding imperfect knowledge. Data on geography, transport, the local and regional political situation, and relative attractiveness of sites can be determined from archaeological survey and historical research, indicating where settlement would have been generally favored and where other regions were likely to be less populated. However, explanations are also needed that demonstrate how settlement transformations, specifically in population, arise, based on uncertainty, while not being overly simplified.
Spatial interaction and structural dynamic models that apply entropy-based (Boltzmann) and dynamic LotkaeVolterra methods (Wilson, 1970;Harris and Wilson, 1978;Wilson, 2012) have the potential to provide explanations that address settlement expansion or contraction within these geographic settings. This paper explores the utility of such modeling and potential insights that might be gained even when information is limited. The goal is to present a generalized simulation that explores how the spatial setting and factors that affect the flow of goods and people can influence urban size transformations and settlement interactions. The paper uses sites surveyed in the KT region in northern Mesopotamia, including those that date to both the MBA and the IA.
We first introduce background information on the data. We discuss both the background of entropy maximizing and structural dynamics modeling approaches. Outputs from the two periods' modeling are then presented. These results explore and demonstrate relevant causal factors that lead to settlement structure differences between the periods discussed. The model and the results' importance to understanding the development of urbanism in the periods are then explored, including broader benefits derived from the approach.

Landscape
The KT is located within the Syrian Jezira, an area measuring an extent of 37,480 KM 2 and confined between the Tigris and Euphrates rivers. The investigated area is largely flat and bounded by the Syrian/Iraqi border to the east, by the Syrian/Turkish border to the north, by the Jebel Sinjar and by the Jebel 'Abd-al-Aziz' to the south, and by the Khabur river to the west (Fig. 1). Scholars have proposed that a major drying trend led to site abandonment by the late third millennium BC (Weiss, 2000), but by the early second millennium BC conditions likely allowed for major settlements to reemerge (Issar and Zohar, 2004).

Middle Bronze Age and Iron Age settlement patterns
Archaeological excavations and surveys carried out in the KT provided the body of data on spatial extent at both regional and local scales as well as occupation histories. While at times such data can be problematic, as sites are obscured from the archaeological record or are simply undetected, considering the scale of the area studied, many sites have been detected because they are mounds that protrude. In the KT, relevant survey data include: Eidem and Warburton (1996), Gavagin (2012), Lyonnet (2000), Meijer (1986), Ristvet (2005), Wright et al. (2007), Ur and Wilkinson (2008), and Ur (2010). Other nearby surveys (Algaze, 1989;Wilkinson and Tucker, 1995;Ball, 2003) have been left out of the analysis, as these are not as continuous with the others.
Within the KT, there are 260 sites that are occupied in the MBA (Fig. 1a). In the western KT (i.e., west of the JaghJagh River), settlements are nucleated and populations were likely concentrated in towns such as Chagar Bazar (ca. 9 ha), Tell Mozan (ca. 35 ha), and Tell Arbid (ca. 7 ha), with the surrounding territory largely devoid of smaller settlements. In the eastern KT, the Tell Leilan (Ristvet, 2005) area alone has 148 sites during the MBA. Here, the dominant role of Tell Leilan is clear, which had an area of ca. 90 ha with many surrounding small villages. Other centers were Tell Farfara (ca. 80 ha) and Tell Muhammed Diyab (ca. 35 ha). Along the Jaghjagh, the main settlements were Tell Brak (ca. 25 ha) and Tell Barri (ca. 9 ha). Overall, far more settlements and greater diversity of site sizes are found in the east; this could be because the area had more favorable climatic conditions (Evans and Smith, 2006). Fig. 1a and Table 1 show and summarize the sites.
In the IA (Fig. 1b), there are 276 settlements in the KT. The west has 194 sites, but only 82 sites are in the east and they were larger (Table 1). In general, settlements during the IA were mostly small and somewhat evenly dispersed in the region (Wilkinson et al., 2005). At the beginning of the IA, few settlements have been located, while in the late IA, during the time of the Assyrian conquest, many small settlements (i.e., less than 5 ha) became established. In the late IA, the aggregate settled area more than doubles from the late second millennium BC (Wilkinson and Barbanes, 2000). During the IA, known larger sites are Tell Halaf (14 ha; Novák, 2009), Tell Hamidya (24 ha;Wäfler, 2003), Tell al 'Id (19 ha), Tell Badan (17.5 ha), Tell Effendi (17 ha), and Tell Beydar (12 ha). The modern city of Nusaybin in Turkey likely overlies a relatively large center, but this site is largely unexplored.
A natural log scale showing settlement size and hierarchies, ranking from largest to smallest, are displayed in Fig. 2. For the western KT (2a), the MBA and IA sites show a significant difference (p-value < 0.5) in site size and hierarchy, using a Kolmogorove Smirnov test, while the eastern KT (2b) there does not seem to be any statistically significant difference (p-value < 0.05) for the two periods' sites. However, this is close to significant (p-value ¼ 0.06).
Overall, comparing all of the MBA and IA sites (2c), the distributions are not statistically significant, but looking at each region and period we see that site size and hierarchy are significantly different (p-value < 0.05). Table 1 shows average and standard deviation in site sizes between the regions and periods. From this, it is clear that the MBA sites are larger in the eastern KT and have greater variation in site size, while in all the KT sites are generally larger in the MBA than in the IA. Table 2 summarizes the top ten largest sties and their size estimates for the periods.

Historical background on the Middle Bronze and Iron Ages
For reasons that are unknown, sometime after 2000 BC, the local kingdoms of Nagar (capital ¼ Tell Brak; Fig. 1 Tunca and Baghdo, 2008) and Mari (Tell Hariri).
During the MBA, the KT is split into various local principalities, not all of which can be localized with confidence. The kingdom of Apûm (Eidem, 2008) seems to control most of the eastern part of the KT, while the rest is divided into several small states including A snakkum, Urke s, and Suduhum (Guichard, 2009) that can be specified with the group designation Ida-maraṣ (Fleming, 2004;Bryce, 2009). For most of that time, these small states accept the sovereignty of the kingdom of Mari (Margueron, 2004), first under the so-called sakkanakku dynasty and then under several rulers that owe their regional importance and their position to their connection with the kingdom of Yamhad (Klengel, 1997), whose centre is Aleppo (Halab) which is outside of the investigated region.
The local states in the KT profit from the important eastewest overland route that leads from Upper Mesopotamia through the KT to the Mediterranean and Anatolia. It is used most prominently by the donkey caravans linking Assur (modern Qala'at Sherqat) on the Tigris in northern Iraq with the region of the Great Salt Lake in Central Anatolia (Barjamovic, 2012). While there is a well established tradition of trade links with the Tigris region, there is a period of about two decades, around 1800 BC, when the KT also comes under the direct political influence of its eastern neighbors. At that time, the Tigris region between Assur and Nineveh (modern Mosul) forms the powerbase of king Samsi-Addu (1833e1776 BC; Charpin and Ziegler, 2003). Tell Leilan comes to serve as Samsi-Addu's local centre, who has the city renamed Subat-Enlil. The dominion of the Tigris region over the KT does not last long after Samsi-Addu's death. For the next decades (Eidem, 2012), the textual sources attest to the military and political maneuvers of regional powers and outside powers from much further away in the  (Charpin, 2012) and Mari is abandoned, the KT seems to become part of the kingdom of Yamhad (Klengel, 1997).
In the early IA, we find the KT at first under local independent rule. Previously, from the 13th to the 11th century BC, the region had been part of the Assyrian state with its power base on the Tigris, centered on the capital city of Assur. But from the reign of A s sur-bel-kala of Assyria (1073e1056 BC) onwards, political control was gradually lost to Aramaean tribal organizations that quickly transformed into regional states and from the 11th century onwards the KT region is fragmented into several small states, with Bit-Bahiani (capital ¼ Tell Halaf; Fig. 1.9; Dornauer, 2010) the most prominent. As the result of successive military offensives during the 10th century BC, the region was once again annexed by the Assyrian state, with its centers further east (Radner, 2011), and the region became integrated into the Assyrian provincial system by 867 BC.
The KT was organized into three provinces. Naṣibina, with the capital city of the same name (modern Nusaybin; Fig. 1), occupied the eastern part of the KT. That area was integrated into the provincial system in 896 BC, first as part of a border march and later as a separate province. This happened probably at the time of the integration of the Aramaean principality of Bit-Bahiani into the Assyrian Empire, which led to the establishment of the province of Guzana by 867 BC. This province took up the western part of the KT and owes its name to its capital city Guzana (¼Tell Halaf). A third province called Raqamatu, after its (unidentified) capital city, which was also known under the Aramaic name Gidara (Bryce, 2009), lies west of Nasibina and north (or perhaps east?) of Guzana. This province was established in 898 BC and may have been merged with another province after 773 BC, when there are no further references to the province (Radner, 2006). Until the fall of the Assyrian state at the end of the 7th century BC, and for more than two centuries, the KT was an integral part of the Assyrian state (Radner, 2002).
What the above historical data indicate is that for much of the MBA we see the KT as politically fragmented and established powers frequently changed, while in much of the IA, particularly when many settlements are well established in the KT, we see a singular large state dominate the region. Below we attempt to explain observed settlement patterns and sizes while reconciling this with the historical data.

Modeling using entropy maximizing and structural dynamics
Entropy maximizing methods have been widely used to model spatial interaction (Wilson, 1970). Here, they are used to represent an indicator of interaction between settlements, which can be considered to represent both migration and trade. Economic or population growth has been modeled using methods analogous to those of Lotka and Volterra in ecology (Harris and Wilson, 1978). The combined models have been labeled as BoltzmanneLotkae Volterra and hence BLV models (Wilson, 2008). In archaeology, various publications have applied similar approaches (Evans and Gould, 1982;Knappett et al., 2008). More recent work on archaeological applications of entropy maximizing methods includes work by Wilson (2012) and Bevan and Wilson (2013). In this approach, we use these and related techniques to model the evolution of settlements in the historical periods in question. Specifically, our aim is to produce simulations of these spatial systems which reflect, as far as we can, the development of these sites, with the goal of explaining why certain sites may have achieved relative prominence and the general settlement size distribution. The validity of the model will be assessed on the basis of the correspondence, in terms of several observed characteristics, between simulated results and data. These characteristics are primarily macro-level features, such as the overall size distribution, so that the agreement is in terms of stylized features rather than highlygranular detail. In the process of seeking such model outputs, examination of those parameter configurations which give rise to the closest correspondence, in this sense, also offers insight into the relative importance of individual factors.
In technical terms, the analysis we present uses a wellestablished formulation of a BLV model, and the primary contribution in this regard is in its adaptation to the present archaeological context. Nevertheless, while the overall structure is standard, a number of methodological innovations have been incorporated which are not present in previous work. In particular, we explore a method for assessing the robustness of results under  uncertainty in data, which is motivated explicitly by the nature of the archaeological data available. The spatial data required for our model is provided in the form of two site datasets; one for each of the periods. For each set e containing 260 and 276 sites respectively e these contain the locations of sites and cost matrices for travel between each pair of sites, calculated via cost surface analysis, using methods applied by Fontenari et al. (2005), in order to account for the effect of varying terrain. We also have rough estimates for the size of each site in hectares, based on surveys referenced, which provide a relative proxy as to which sites likely had larger populations.

Model formulation
The model takes the form of a standard spatial interaction model (SIM) of the type which has been used previously in other contexts (e.g., Harris and Wilson, 1978;Bevan and Wilson, 2013), with some modifications for the present setting. We define the following variables for each of the sites: , people and goods) originating at a given site i Z j ¼ 'attractiveness' or size of site j; a variable that represents any factor that makes a site attractive to settle d ij ¼ distance (i.e., cost of travel) between any two sites i and j, normalized by the mean of all such distances r ij ¼ number of rivers crossed by a direct path from i to j In the case of cost values d ij for which i ¼ j (representing, notionally, the distance from a site i to itself), these are set to f, a parameter which represents impedance to internal trade at a given site, before the normalization step. In addition, for every river which is crossed by a direct path between two sites i and j, an additional impedance r is added to the corresponding cost value, so that the total is (d ij þ rr ij ).
One further factor worthy of note here is that of external interaction; that is, each site's interaction with sites outside the area of study. This would be manifested as an additional in-flow for each site and therefore contribute to site growth. A simple representation of this in the model would be the inclusion of a constant term for 'rate of growth due to external interaction', encapsulating all such effects. Simulations were carried out with such a term included, but it was found that its effect could also be achieved, in general, through variation of f. This can be rationalized by considering the function of f as a barrier to internal interaction: such interaction is a feedback effect, providing in-flow to a site in proportion to its current size, and so increasing this flow mimics the consequence of external influence.
The model itself can, broadly, be divided into two parts: the estimation of interaction flows and the use of these flows to determine site size evolution. The first of these relies upon establishing, for any pair of sites i and j, the utility of interaction with j as perceived by i. This takes the form of a cost/benefit calculation, with the benefit being a function of the size of j and the cost representing the physical impedance (i.e., distance between). Two parameters, a and b, are introduced at this stage: a specifies the scaling of utility as size varies, and b determines the strength of the negative effect of cost. Using entropy-maximizing methods, the most likely set of flows, given that the total flow originating at each site is known, is then found. The allocation given by this method is relatively simple: for any given site i, the total out-flow X i is shared between each other site j in proportion to the utility of j. In this way, we find the flow S ij between any such pair of sites.
In the second part of the model, these flows are used to update the site sizes; that is, to determine the growth/decline of each site under the calculated set of flows. For each site j, the total inward flow is found by summing S ij over all i, and this value is compared with the site's current size Z j . If the total flow is less than Z j , the current size is unsustainable and the site shrinks, whereas a surplus of inward flow leads to the growth of j. The model, in essence, is appropriate for measuring generic factors (e.g., environmental constraints, political conditions, economic conditions) that could lead to some sites becoming larger than others.
Having made these calculations and adjusted X i and Z j accordingly, the model proceeds to the next time-step and the process is repeated. For any individual simulation, we run the model for a set number of discrete time-steps dt: The flow S ij between each pair of nodes i and j is calculated using the following formula, which utilizes a and b: These flows are summed to give the total incoming flow D j to each site j: This incoming flow is used to calculate Z j at the next time step, Z j (t þ dt) , via the formula: We now find X i (t þ dt) for the next time step by taking it to be the corresponding Z i (t þ dt) value, normalized according to the total of all Z i (t þ dt) and rescaled so that the X i (t þ dt) continue to have mean 1: Then return to (1).

Model outputs
Having completed any given numerical simulation of the model, a natural question arises of which quantities should be regarded as the fundamental outputs of the model; these are precisely those quantities which will be used to evaluate settlement sizes and distributions for the periods discussed. For each simulation, we have several observable features upon which to base assessments: the final site size Z i for each site j, and the distribution of these values across all sites. Since we have numerical estimates for these values within our datasets (the distributions of which are shown in Fig. 3) we can make meaningful numerical comparison between these distributions using statistical tests; the structure and properties of a NystueneDacey (NeD) graph (Nystuen and Dacey, 1961; construction described in Appendix A) derived from the simulated flow matrix S ij . Such a graph defines a hierarchy of sorts for sites, with links representing subordinate relationships. This facilitates the analysis of both graph-theoretical properties and the average geographical length of the estimated links, i.e. the typical length traveled to a locally-dominant site. Although true flow data is not available for comparison, we can compare the link length with that of another graph derived from the real data, for example by linking each site to its closest 'large' site; The sites are colored according to size, ranging from yellow to red as size increases, and the links are those of the induced NeD graph. It is notable that the links are significantly longer in the low-b simulations where there is a lower associated cost to travel, and that the changing value of a has the effect of picking out different dominant sites. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

T. Davies et al. / Journal of Archaeological Science 43 (2014) 141e154
the extent to which the largest sites in the real data are found to be large during simulated runs. 'Large' is taken here to refer to sites whose size is an order of magnitude higher (i.e., 10 times) than other sites: for each dataset, we identify a set of L such sites (L is 19 and 6 for MBA and IA respectively) and, for each simulation, record the number of these which appear in the top L simulated sites, when ranked; visual inspection of outputs and evaluation of their historical validity.
It is in these terms, therefore, that we characterize individual simulations; however, we stop short of identifying specific criteria for what constitutes a 'match' with real data. Firstly, each of the above points represents a well-justified objective for a model and it is far from clear how they should be combined. More importantly, though, to be prescriptive in this way implies a level of confidence in the accuracy and completeness of the data, which is not present, and therefore to fit in this way may be counter-productive.

Results
There is clearly a wide degree of variation possible within the model, in terms of parameter choices, manipulation of the underlying datasets, and synthesis of results. In the following sections we explore several of these, incrementally building the complexity of the approach. All simulations are run for 3000 ticks, with step size dt as 0.01; in all cases this was sufficient for the site sizes to reach a steady state. Initially, all variables were the same for sites, while only distances (i.e., cost surface) between sites were different. Data, code, and instructions for the model are made publically available (http://discovery.ucl.ac.uk/1414597/).

Parameter variations
The model features four parameters e a, b, r and f e which can be varied to give different outputs, for a given set of initial conditions. Each has a real-world interpretation: a represents returns to scale of attractiveness as a function of settlement size, with a value of 1 indicating a linear relationship; b determines the impedance of distance (i.e., how easily goods and people can move across given distances); r governs the extent to which rivers inhibit interaction; f controls the relative attractiveness of internal and exogenous activity for a site, which accounts for edge effects as this addresses outside interactions and sites not addressed in the model.
The effect of such parameter variation can be seen visually in Fig. 3, which shows outputs for the MBA as a and b are varied. In these simulations, the system is initialized in a homogeneous state; specifically, X i (0) is set to 1 for all sites i. These particular simulations are intended not as an attempt to produce realistic results, but instead to give a visual impression of the effect of parameter variation: it is clear to see that b is a significant determinant of the length-scale of interactions, while a influences which sites become dominant.
A different perspective can be seen in Fig. 4, which shows the average length of connections in the induced NeD graph for the IA, across a range of b and r. While this offers limited insight into the underlying system, these results exemplify the significant and highly non-linear manner in which the model parameters interact.
As would be expected, smaller b values e representing lower impedance to travel e lead to longer interaction links. The behavior with r, however, is more interesting, particularly at low values of b: as impedance increases, link length decreases, before unexpectedly increasing towards higher values. An explanation might be found by considering the relative positions of rivers and sites: in the southern area, many sites are clustered around rivers, and the impedance of the river increasing above a certain value might cause their dominant interactions to switch from those involving river crossings to favor sites on the same side of the river, which may be further away.

Fitting to empirical data
The ultimate objective of analyzing parameter space in this way is to identify a set of values that correspond most closely to the real data. Although this correspondence could be defined in many ways, we take as our first objective the distribution of site sizes. This requires a measure of the similarity of distributions, and for this we use KolmogoroveSmirnov (KeS) distance. For any set of data, it is possible to calculate an empirical cumulative distribution function (ECDF): a step-function which, for any value x, gives the probability of finding a value less than or equal to x when randomly sampling from the data. The KeS distance between any two distributions is simply the maximum difference between their ECDFs, over all possible x; it therefore ranges between 0, in case of a perfect fit, and 1 in the worst case. For each output, we compare the distribution of the modeled site sizes with that of empirical site size estimates, having first normalized each set by dividing by its mean.
In order to investigate the goodness of fit for various configurations, a grid search of parameter space was carried out and the Ke S distance found in each case. One feature which becomes immediately apparent is the disparity between the two periods in terms of the parameter configurations which give rise to the closest correspondence. Fig. 5 shows the results of simulations for a sampling of the parameter space of a, b, r and f, shown as an interpolated 'heat map.' Presented in this way, the difference in character between the two periods is immediately apparent, since the areas of parameter space indicating good agreement differ markedly between the two plots. For the MBA, better fit is found for higher values of a, for example, implying that returns to scale were larger for this system. In addition, higher values of f also lead to closer fits in MBA than in IA, suggesting that the MBA sites were less reliant on internal interactions. Above all, the dissimilarity between the plots suggests a significant difference in character between the two periods.
Using the results of this parameter search, configurations were found for MBA and IA sites, which give the minimal KeS distance. The relevant parameters are given in Table 3, and the results are shown in Figs. 6 and 7, which show site size distributions and mapped results respectively. The relatively good fit e seen both visually and in the low KeS values e indicates that, on this gauge, the model has successfully reproduced the empirical distribution.
Although these results show good agreement in terms of site sizes, the sites which have large estimated size in the observed data tend not to be those which the model forecasts to be large. The failure to identify precise sites does not necessarily imply that the model has not identified more general geographical areas, which would be expected to feature large sites. Of course, in reality, there are various factors that lead to the dominance of certain sites (e.g., earlier settlements could initially be larger, socio-environment factors may favor one site vs. others), whereas only geographic location is included here. Considering the case of Tell Leilan, for example, it is clear from inspection that there is no a priori reason why it, rather than its nearest neighbor Rehaya, should become dominant in the model. If the analysis is widened, though, to consider situations where the dominant simulated site is not one of the true largest sites, but is in close proximity to one, it is possible to find the extent to which the localities of large sites are identified.
To do this for any given simulation output, we take each of the 10 largest simulated sites and calculate the distance to the nearest of the 'true' 10 largest sites (that is, the minimum distance one would have to travel to reach one of the true top 10). The average of these values is an estimate of the extent to which the true large sites have been 'missed', in geographical terms. In the MBA case, this value, for the parameters given in Table 3, is 7564 m. To put this value in perspective, we can compare it with the value it would take if, instead of using the top 10 simulated sites, a random set of 10 sites were chosen and the distances from these to the true top 10 calculated. Doing this 100 times, we find a mean value of 15,597 m with standard deviation of 4648; that is, the sites which are forecasted as large by our model are significantly closer to the true large sites than would be expected by chance. Although comparison with random points is relatively unambitious, this value, along with visual comparison, suggests that the model is indeed finding sites that are geographically close to known large sites, suggesting the model is far better at forecasting the largest sites than random selection. The equivalent value for the IA is 4272 m, with a corresponding null value of 12,062 m (s.d. 2602).

Initial conditions
The results of the previous section suggest that, when we seek to produce a distribution of site sizes similar to the real data, known large sites do not, in general, emerge in the model. In order to distinguish such sites, therefore, some modification to the model is required. A natural approach to this, within the current framework, is to manipulate the initial conditions for site size; rather than taking these to be uniform, we adjust these according to what is known about the sites in question. Such an adjustment perturbs the model in favor of certain sites; in more concrete terms, this could be regarded as accounting for any effects not explicitly included in the model which currently only considers the spatial configuration of the system.
Our approach to this is to divide the sites into three categories, according to the order of magnitude of their given size in the data (essentially a classification as 'small', 'medium' or 'large' on the basis of their hectare sizes lying in the ranges [0,1), [1,10) or [10, N), respectively). The initial condition for each site is then set to be the mean of the known size of all sites in its category. In the MBA, for example, the mean size of 'large' sites is 28.73, and so all large sites are initialized with this value. Although guiding the system, this still represents a relatively coarse-grained classification of sites. Fig. 8 shows the result of simulations under this modification, with parameters identical to those shown in Table 3. In terms of top sites, a modest improvement is apparent: for the MBA, the sites at Hansa, al-Andalus and Khirbet el-'Abd are now significant, and in the IA there is now a larger site at Tell Halaf.
If we allow ourselves to make changes to the parameters, however, the results can be improved still further in these terms.
For the MBA, using the values a ¼ 0.7, b ¼ 45, r ¼ 1000 and f ¼ 500 leads to the output shown in Fig. 9a. Several empirically large sites, such as Tell Brak, Tell Hamidiyah and Tell Leilan, are now amongst the largest forecasted, and several of the other sites found to be large (e.g., Duger, Lazzaga and Khirbet el-'Abd) are known to be significant. For the IA, on the other hand, the values used are a ¼ 0.4, b ¼ 45, r ¼ 5000 and f ¼ 500, with the output shown in Fig. 9b. We can, of course, also consider the distributions of site sizes in these cases: we have KeS values of 0.21 and 0.19 respectively. Although these do not represent such a close fit as the previous section, the difference is not so large as to nullify the results.
The ability of the model to identify these sites is encouraging: it offers further confirmation that the model is capable of producing sensible outcomes and demonstrates that only a small amount of additional information is required to improve the accuracy of the simulation. It is, however, instructive to consider the nature of the changes to the model that were required to give rise to such results. The need to change initial conditions suggests the existence of certain site-specific factors, which have a material effect on the real-world outcome but are not included in the basic model (which incorporates geographical location only). Without accounting for these factors in a more elaborate model, the results of this section, achieved by imposing the three-tiered initial classification, should be viewed as illustrative only.
In terms of parameters, we may draw some insight from an inspection of results e suggesting that the most crucial change in terms of the identification of major sites is the lowering of f. This can be rationalized by considering the spatial layout of sites in both periods, which includes significant clustering in both cases. Since every site represents a source of flow, a site which is situated in an area of high site density has a large potential in-flow within relatively short distance. If such a site is able to dominate locally, as tends to be the case within the model, it is therefore likely to reach larger size than a similarly dominant site in a sparse region, simply because of the greater resource on which it may draw. Lower values of f represent greater emphasis on internal development at a site; as noted previously, this might also be thought of as giving greater weight to extra-regional trade and, therefore, mediate against the above effect. A good example of this is the southern region around Uwayja in the IA period, which has a high density of sites: in Fig. 8b, this density is clearly associated with large site sizes, whereas the effect is much less pronounced in Fig. 9b.

Partial dataset sampling
The clustering of sites mentioned above may be the result of unevenness in the surveying of sites, but also relates to a more general problem with the data: that of the contemporaneousness of sites. The possible lifetimes of the locations we are modeling lie within a wide range; at any given time, it may be the case that only a subset would be in existence (and these subsets may not be known). To ameliorate this, we use a new simulation protocol, which involves repeated sampling of the sites.
To do this, instead of performing one simulation run for a given parameter configuration, we perform a series of M such simulations. For each simulation, we randomly sample a subset of R of the sites,    Table 3, but with initial conditions modified according to the known size of sites, as described in the text.

T. Davies et al. / Journal of Archaeological Science 43 (2014) 141e154
and carry out the simulation using only these sites. When all simulations are complete, we take the mean modeled size of each site across all the simulations in which it featured, and take this to be our estimates of that site's size. The underlying principle of this method is that each sampling represents a possible 'state of the world', and that averaging in this way accounts for our uncertainty about the composition of the dataset: when a site is consistently modeled to be large, for example, we can be confident that it is likely to be large whatever the true composition of the system should be.
We carry out this process for each of the MBA and IA datasets, taking M ¼ 500 and R ¼ 100 and with the other parameters as specified in Table 3. Encouragingly, there is good consistency between the results of these simulations and those of the original simulations, which used the whole dataset. We can use Spearman's rank correlation coefficient (which gives a value of 1 in the case of perfect correlation and 0 in case of total absence) to examine how closely the rankings of the sites under each protocol match; this quantity is 0.556 for the MBA and 0.976 for the IA. This is significant in both cases, and indeed remarkably high in the IA case, suggesting consistency of results. In other words, sites which are found to be large under the original algorithm remain large under this averaging and random sampling system, and vice versa.
If we consider the variation of individual sites across such a run of simulations, we observe a further difference between MBA and IA periods. Fig.10 shows the sizes of individual sites across various runs; it is noticeable that the distributions are normal-like in the IA case, whereas the MBA case has a much more flat distribution, with many low values. This can again be understood in terms of the clustering of sites: many IA sites are well-clustered and are therefore somewhat interchangeable from the point of the model, whereas greater spatial diversity of the MBA sites means that system is much less resilient to changes. The results are also shown in map form in Fig. 11.
In terms of size distribution, however, those simulations show a significant disparity with the reference distribution (in terms of Ke S distance); in particular the distribution has much lower variance. This can be understood in terms of the averaging procedure; any given site will be modeled as large in some of the simulations and smaller in others (as shown in Fig. 10). Averaging over runs gives an impression of which of these regimes it tends to be in, but does not necessarily reflect its true size. Under this principle, the true size is the outcome of the model under one particular (but unknown) sampling; the averaging simply offers insight into its tendency across all possibilities.

Discussion
The results demonstrate the utility of the model as a method that explores spatial systems in historical contexts such as those considered here, while also forecasting, at least far better than random sampling, lager sites nearer to their true location using input survey data. The model produces important insights and demonstrates broader utility in describing settlement hierarchies and detecting a variety of factors that may have caused some sites to become relatively larger than others. Results show potential political and historical differences affecting settlement in the MBA and IA that can be explained by the advanced model.
Several aspects are particularly notable, both from the point of view of modeling and of providing insight into the underlying system. Firstly, the model is able to produce close matches with real data in terms of the distribution of site sizes, even in the most basic implementation that focuses on spatial location. As well as providing partial validation of the model, this stage of the work is notable for the varying parameter configurations that are required to produce the closest matches in each time period. The fact that the best fits are found in markedly different areas of parameter Fig. 9. Mapped simulation outputs using 'small/medium/large' initial conditions as described in the text, with parameters chosen so as to identify key sites: a) MBA with parameters space e particularly in terms of a and b e implies that there is an inherent difference between the periods. The results indicate that returns to scale in site size were greater in the MBA (due to higher a) and that interactions typically had a shorter length-scale (higher b). In the IA, lower a indicates less importance on site size, while lower b indicates greater movement or facility of movement across the landscape. Translated to the discussed historical data, a politically fragmented landscape, such as in the MBA, could restrict movement, thereby leading to larger local sites that also absorb more flow from nearby sites. There may have been more emphasis on establishing large sites, perhaps for protection (i.e., which reflects higher a) during periods of political competition, further encouraging larger sites to become even larger.
We can interpret the results by the fact that in the MBA the KT is characterized by the presence of several peer polities aiming to exert their influence over their surrounding hinterland, while in the IA settlements are well integrated into the Assyrian imperial structure. During the MBA, we see the Jaghjagh River playing an important physical/political boundary, and modeling results demonstrate how river impedance can differentially affect site size distributions depending on which side of the river bank a site might be located. While in the IA, the landscape would have been easier to navigate due to the Neo-Assyrian empire's control of the region, with movement likely to be less restricted, particularly for a long period of time (i.e., more than 200 years). Easier movement enables the dispersion of flow to many sites, leading to greater similarity in site sizes. Diminished a indicates less of a need or ability to have larger sites in the IA, which could be caused deliberately by the Neo-Assyrians or other environmental factors.
Overall, and another important result demonstrated in the model, is that it was far better at forecasting larger sites near their true locations than by selecting sites randomly. In order to locate such sites, it was necessary to manipulate both the initial conditions of the model and some of the parameter choices. The first of these implies that it was necessary to include a certain amount of known information about the sites in question; in other words, that location data alone were not sufficient to explain their dominance. Nevertheless, the manipulation of initial conditions here was done in a coarsegrained way, suggesting perhaps that any extra information need not be particularly detailed. In the case of parameter variation, it was found that closer matches tended to be associated with lower values of f, a parameter controlling local and exogenous interactions, which can be related specifically to the nature of the data. The fact that the data contained several areas of high site density presents something of a modeling challenge. A large number of sites translates to a high activity potential in the model, which may not represent reality: the high density of sites may simply be an artifact of the area being relatively well-surveyed or other distortion in identification of sites.
The parameter f has the potential to reduce this effect, by placing greater interpretation on internal site features or interaction with exogenous, non-KT sites outside of the simulated area, depending on interpretation. In a sense, then, it is able to 'tune' the relative influence of local interaction effects and address edge effects.
The results generally show that other non-geographic factors do make specific known sites larger, as sites become relatively large through modifications to initial size and f. This can, for instance, reflect external political conditions or initial site advantages. In particular, it seems that Tell Brak, Tell Leilan and Tell Mozan in the MBA would need greater exogenous effects or initial advantages to enable them to consistently reach a relatively greater size. This could come in the form of Samsi-Addu making Tell Leilan (the ancient Subat-Enlil) an important political capital, Tell Brak being the seat of the "Lady of Nagar" (Oates et al., 1997, 141), and local or external benefits that Tell Mozan may have relative to other sites. In the IA, the Neo-Assyrians choosing Tell Halaf as a provincial capital could enable that site to be relatively large.
The issue of clustering of sites is once again apparent when we consider another method of simulation: that of running repeated simulations on only partial datasets (i.e., see section 4.4). The motivation for this is to combat the effect of temporal uncertainty in the data, since possibly only some of the sites are likely to be contemporary. The IA case, remarkably, shows little difference when using this method, primarily because of the spatial structure: high clustering means that the role of deleted sites is essentially assumed by close neighbors and in simulation runs site size hierarchies are very comparable. This indicates that even if some sites were occupied at different times or missed in survey, the hierarchy of site sizes would not be significantly different than that achieved in the results. The situation differs in the MBA, which appears to be more volatile in this sense and consistency of larger sites and hierarchies are not as easily maintained across various runs, although the results are still fairly consistent based on the Spearman's rank correlation coefficient. Nevertheless, based on this uncertainty, this serves to emphasize that small changes in this dataset are liable to have a significant effect on model outputs.
Despite the promise of the results, some shortcomings of the modeling approach must be acknowledged. While also being a strength to the modeling approach, as it more easily allows the model to be transferable to other case studies, the model is necessarily reductive and implements several factors in a relatively coarse manner. The parameter f, which encapsulates a number of factors, is a case in point, and is a consequence of seeking a parsimonious model. In a similar vein, few non-spatial factors are included, and these may have significant influence on the true evolution of sites; the failure of the model to pick out certain sites can be ascribed to this. It must be emphasized that the fine-grained accuracy of the model is directly related to its complexity, and that the work presented here is a test of the potential of a simple model on imperfect data. In future developments, several of these issues could be addressed. In addition, the question of the selection of optimal parameters remains a subject of ongoing research in light of the incomplete nature of archaeological data.

Conclusion
The goal of this paper has been to present a method that explores how the spatial setting and other unspecified socioenvironmental factors in the MBA and IA could affect the flow of goods and people that ultimately influence urban growth and settlement sizes. We explore certain aspects from the survey data and try to match known historical events to explain modeling results. The strength of the modeling approach is that it can explain the role of such factors as location, settlement importance, exogenous effects, and movement and flow of goods and people on site sizes and their distribution. While we are not completely certain what caused observed results, the model is reconcilable with known political and historical events. In terms of how outputs might be used, there are several options for how results might be used for meaningful historical/archaeological understanding. On the one hand, we can use point estimates, using optimized parameters as we have done here, to forecast the relative size of known (or unknown) sites. It is clear from our results, however, that there is a range of parameters which all lead to very similar outputs, and that selecting only one of these might not be prudent. An alternative to be explored would involve the establishment of some criteria for the 'acceptability' of a simulation output, identification of all parameter configurations that satisfied these criteria, and the production of an output which represents an average across these configurations.
The present work has potential to be extended in various directions. On the one hand, there are various aspects of the model where it should be possible to incorporate uncertainty to a greater degree. For example, our consideration of known site sizes uses point estimates for these, where it may be more prudent to draw values from a weighted distribution. In the case of rivers, also, the impedance which we associate with river crossings, r, is constant everywhere, which represents a simplification. Given more detailed data, it would be possible to experiment with varying this value on a river-by-river basis, or by incorporating known crossings. On a deeper level, given that the clustering of sites in the data has been identified as a significant complicating factor, further work is likely to be required to consider how this problem might be tackled. This might be done by developing a process whereby sites Fig. 11. Mapped results for a) MBA and b) IA under the repeated sampling procedure. Sizes shown here are the average sizes for each site across 500 simulations, each of which includes 100 randomly selected sites only. The parameters used are those given in Table 3. which are consistently insignificant can be deleted, or merged with other sites, in order to reduce imbalances of density. As for larger benefits to archaeology, the model demonstrates that it can begin to reconcile socio-political or even environmental constraints to an abstract level, whereby such circumstances can be translated to parameters such as a, b, and f that affect the importance of site size, cost of movement, and internal-exogenous site interactions respectively. While such variables can be affected by a host of circumstances, the benefit is the method can be transferred to other and very different case studies more easily. The model is applicable to describing observed settlement structures by providing parameter estimates that enable nearly matching simulated settlement data. The model is also useful for forecasting general areas where larger or smaller sites are to be expected, demonstrating its potential for helping archaeologists to locate such sites. Overall, such modeling allows archaeologists to differentiate how significantly different periods are in site distribution and size structures and begin to understand this with causal factors.