Integrating Global Sensitivity Approaches to Deconstruct Spatial and Temporal Sensitivities of Complex Spatial Agent-Based Models

: Spatialagent-basedmodels(ABMs)canbepowerfultoolsforunderstandingindividualleveldecision-making. However, in an attempt to represent realistic decision-making processes, spatial ABMs often become extremely complex, making it difficult to identify and quantify sources of model sensitivity. This paper im-plements a coastal version of the economic agent-based urban growth model, CHALMS, to investigate both space- and time-varying sensitivities of simulated coastal development dynamics. We review the current state of spatially- and temporally-explicit global sensitivity analyses (GSA) for environmental modeling in general, and build on the innovative but nascent efforts to implement these approaches with complex spatial ABMs. Combined variance- and density-based approaches to GSA were used to investigate the partitioning, magni-tude, and directionality of model output variance. Time-varying GSA revealed sensitivity of multiple outputs to storm frequency and cyclical patterns of sensitivity for other input parameters. Spatially-explicit GSA showed diverging sensitivities at landscape versus (smaller-scale) zonal levels, reflecting trade-offs in residential housing consumer location decisions and spatial ‘spill-over’ interactions. More broadly, when transitioning from a conceptual to empirically parameterized model, sensitivity analysis is a helpful step to prioritize parameters for data collection, particularly when data collection is costly. These findings illustrate unique challenges of and need to perform comprehensive sensitivity analysis with dynamic, spatial ABMs.


Introduction
. Spatial agent-based models (ABMs) are capable of investigating the emergence of complex phenomena, such as land-use change, that result from many distributed interactions between heterogeneous agents and their environment, which are not well represented by top-down statistical or general equilibrium modeling approaches (Irwin et al. ; Huang et al. ). Many spatial ABMs are characterized by multi-level and nonlinear interactions, emergent behavior, and path-dependent outcomes (Brown et al. ; Ten Broeke et al. ). While these features are o en desirable to capture salient features of real systems, model complexity makes it di icult to deconstruct sources of model output sensitivity in relation to variation in inputs. This task is particularly important when transitioning from a conceptual to empirically parameterized model. When data collection is costly, such as attaining behavioral data for ABMs from surveys, sensitivity analysis is a beneficial step to identify parameters responsible for model sensitivity and prioritize data collection e orts to constrain those variables. This paper demonstrates the synergies of integrating two di erent methods for sensitivity analysis to achieve a comprehensive investigation of spatial and temporal sensitivities of an established conceptual ABM for moving towards a more empirically grounded model.

.
Standardized methods for sensitivity analysis for ABMs are still developing (see Ten Broeke et al. ; Thiele et al. ), but there is a wealth of knowledge to draw from in the broader environmental modeling field. Global sensitivity analysis (GSA) is a relatively well-developed methodological area for investigating how variation in model inputs leads to model output variation (Brugnach ; Cariboni et al. ; Marino et al. ; Pianosi et al. ; Sarrazin et al. ). GSAs are implemented through a set of mathematical techniques designed to quantify the relative contributions of and propagation of uncertainty from individual model parameters and their interactions with other parameters through a comprehensive sweep over the parameter space (Ligmann-Zielinska & Sun ; Pianosi & Wagener ). This is distinct from simple one-factor-at-a-time (OAT) methods, which are more suitable for exploring individual parameter e ects on outcomes (e.g., scenario analysis) (Ten Broeke et al. ). A wide array of GSA methods exist that range from simple parameter perturbation to variance decomposition, and are designed for equally diverse purposes ranging from model fitting to understanding the robustness of model outcomes across parameter values (Pianosi et al. ; Ten Broeke et al. ). GSA methods have recently been reviewed in detail in by Pianosi et al. ( ). A variety of GSA techniques for specifically analyzing spatial and/or temporal model sensitivities have also been developed and reviewed elsewhere in detail (Saltelli et al. , ; Lilburne & Tarantola ). However, these GSA methods may not be su icient to address the unique challenges of ABMs (Ten Broeke et al. ), and the few recent e orts to apply GSA methods to spatial ABMs demonstrate the inadequacy of any single existing GSA method (e.g., Lilburne & Tarantola ; Ligmann-Zielinska & Sun ; Saltelli et al. ; Ligmann-Zielinska et al. b,a; Thiele et al. ). Thus, understanding of the sensitivities of many spatial ABMs in use today remains limited (Ligmann-Zielinska & Sun ).

Variance-and density-based GSA .
Two broad approaches to GSA, variance-and density-based, have been most o en applied to spatial environmental models. Variance-based GSA quantifies the sensitivity of outcome variance attributed to variability in model inputs, and has the advantage of being both model independent and relatively easy to implement through the calculation of variance sensitivity indices (Ligmann-Zielinska & Sun ; Saltelli et al. ; Pianosi & Wagener ). However, the generation of su icient model output to obtain accurate sensitivity estimates with variance-based approaches can be computationally expensive, particularly for sophisticated models with high levels of parameter interactions (e.g., Baroni & Tarantola ) required , model evaluations for uncertain inputs. Additionally, variance-based analysis assumes output variance is a reliable sensitivity measure, but this may not be the case when model outcomes display skewed or multi-modal distributions (Pianosi & Wagener ).
. Non-normal outcome distributions are commonly observed in complex systems due to a wide range of nonlinear, spatial processes (e.g., agglomeration), such as economies of scale or path-dependence of development location due to infrastructure provision (Brown et al. ; Manson ). In these situations, variance-based GSA may be misleading, but density-based methods for GSA, which are 'moment-independent', are useful. Model sensitivity is determined by comparing model output distributions generated by di erent input parameter values, rather than comparing specific moments of output distributions (i.e., variance) when variation in one input is fixed (Pianosi & Wagener ). Pianosi & Wagener ( ) have developed a density-based sensitivity index, PAWN, which uses di erences in cumulative distribution functions (CDFs) of model output with both variable and fixed input values to quantify model sensitivities. A particular advantage of this approach is that model sensitivities can be investigated within particular parts of a given model output distribution, which is important for contexts with low probability-high impact events. While density-based approaches cannot (as of yet) provide insight into the dynamics of model sensitivities, they can validate and complement variance-based approaches in situations with model output that is not normally distributed (Pianosi & Wagener ). However, density-based approaches have seen limited use in the environmental modeling domain and have not been applied in a spatial context.

Space-and time-varying GSA .
The GSA methods reviewed thus far are most commonly used as 'final snapshot' analyses of model outcomes, which are inadequate to fully explore sensitivities of complex ABMs (Richiardi et al. ; Ligmann-Zielinska & Sun ). Analyzing only the final state of model outcomes does not provide insight into the relative and varying importance of particular parameters throughout model execution, and potential for synergistic interactions between model parameters leading to locally or globally nonlinear behavior (Ligmann-Zielinska & Sun ). Furthermore, temporally and/or spatially distributed inputs and outputs make GSA based on Monte Carlo methods used for conventional environmental models are conceptually challenging and computationally expensive to implement for ABMs (Lilburne & Tarantola ).
. Spatial GSA has received considerably more attention than time-varying GSA, and has been led by the geographic information sciences (GIScience) and hydrological modeling communities. Early e orts to analyze spatial sensitivities of models relied on methods of simplifying space in order to make analyses more computationally tractable. For example, Hall et al. ( ) divided the model's spatial domain into distinct zones, each characterized by spatially aggregated scalar inputs, to which GSA techniques were independently applied. In cases where input parameters vary across space, GSA may be applied to each individual grid cell or point (e.g. Tang et al. ), but the e ects of spatial interactions on model outputs may not be captured (Lilburne & Tarantola ). An alternative approach entails randomly or purposely selecting a subset of grid cells or points on which GSA can be carried-out (e.g. Avissar ; Dubus & Brown ; Lilburne et al. ). In this sub-sampling approach, a cluster analysis can be used to identify spatial regions that are internally consistent and distinct from other regions and may therefore display di erent sensitivities (Lilburne & Tarantola ).
. Progress in time-and space-varying GSA for ABMs has been advanced largely by Ligmann-Zielinska and colleagues. Ligmann-Zielinska et al. ( b) introduced a multi-scale, variance-based approach to GSA in which sensitivity analyses were conducted to understand variability in modeled nutrient loadings at the regional level as well as for individual lakes. The regional scale GSA identified the input factors, such as fraction of agricultural land converted to fallow and farmers' decision-making rules, which where the sources of regional variation, while the lake-level analyses demonstrated sensitivity to watershed-specific parameters, such as value of production and prioritization of land characteristics. Ligmann-Zielinska & Sun ( ) examined time-varying GSA of a land-use ABM to understand how the contribution of particular input to model output uncertainty and sensitivity waxed and waned through model execution. They found that spatial outcomes, such as the size of the largest contiguously developed area, showed varying sensitivity throughout model execution due to feedbacks among landscape characteristics, such as perceived scenic beauty and land values, that varied as development patterns changed.

.
This paper advances previous work by integrating time-and space-varying GSA using variance-and densitybased approaches in tandem. We illustrate this integrated approach for ABMs through application to a coastal version of a generalized economic urban growth ABM, CHALMS (Magliocca et al. , ). Coastal land-use patterns are particularly complex due to interactions between consumer locational preferences and risks of storm-induced property damages, which are influenced by processes that are o en disconnected in time and space. ABMs have the capability to provide insights into the trade-o s between coastal amenities and storm damage risks individuals might make in their location decisions (e.g., Filatova et al. ; Filatova & Bin ), but parameterizing the decision rules (e.g., utility function) in such context is di icult and can be a main source of model uncertainty. Thus, we use a suite of GSA methods to investigate whether sensitivities are more pronounced at particular times in a coastal region's development and/or in particular locations on the landscape (e.g., along the coast versus further inland). These analyses determine the extent to which input parameters for coastal amenity values, homeowner location preferences, and storm frequency and severity settings contribute to variability in the area developed, number of houses built, and diversity of the housing stock over time and across several 'spatial zones' within the landscape. As with all GSA methods, results do not establish causal relationships between model inputs and outputs. In the case of variance-based GSA, directionality of model output variance cannot be established. For density-based GSA, information about dynamics of model output variance are lost. However, combining the two types of GSA approaches can identify key parameters responsible for sensitivity of model outputs and, in some cases, indicate the direction, timing, and spatial structure of model sensitivities in response to parameter perturbations -thus acting as a proxy to study model causality. Such outcomes are essential for advancing model development from a conceptual to empirical spatial ABM. .
The remainder of this paper proceeds as follows. The next section provides a brief overview of the coastal CHALMS model (Section ). The following sections describe the time-and spatially-varying GSAs that are carriedout by varying seven key model inputs (Section ). Spatial and temporal model sensitivities are then presented, followed by a discussion of the benefits of and limitations to GSA for spatial ABMs.

Coastal CHALMS overview
. The ABM presented here is a version of the CHALMS model (Magliocca et al. , ) adapted to a coastal landscape. The model simulates the conversion of undeveloped land to residential housing through spatial and utility-or profit-maximizing decisions by residential housing consumers, a representative developer, and landowners ( Figure ). A full model description using the Overview, Design concepts, and Details + Decision-Making (ODD+D; Grimm et al. ; Müller et al. ) is provided as a supplemental document. The main objective of the original version of CHALMS is to understand which landscape features and agent interactions are most important for explaining the timing and location of residential housing development. In the coastal version presented here, the influence of coastal amenities and risk of property damage from storms is also considered. To do this, CHALMS explicitly links housing and land markets to simulate endogenous, coupled market prices. Understanding spatially explicit and time-varying dynamics that result from land and housing market interactions is essential for evaluating potential outcomes of land-use policies (e.g., zoning; Magliocca et al. ).

Landscape and model initialization .
The simulated landscape is a stylized representation of a developing coastal region, i.e., does not represent any particular geographic area. Land acquisition and, in part, residential location decisions are based on two landscape attributes: distance to a central business district (CBD), which influences travel costs to places of employment, and proximity to the coastline from which consumers derive utility in the form of a natural amenity ( Figure ). Potential damages from storm impacts are also related to proximity to the coast (decreasing with distance from the coast) and the annual probability of storm occurrence. Travel costs increase linearly with distance from the CBD at a rate of $ . per cell ( ) decay nonlinearly as distance from the coast increases ( Figure ). Spatially-varying expected costs associated with property damages from coastal storms are a function of both proximity to the coast ( Figure ) and the annual probability of storm occurrence (Costanza et al. ). The model is executed for discrete time steps (one model time step equals a year) over a gridded cellular landscape of , cells ( by ) and each cell is equivalent to one acre. The simulation time span was chosen based on the assumption that the development process would experience minimal fundamental changes over years (e.g., from financing innovations), and the spatial extent was chosen to minimize any spatial boundary e ects on development given the simulation time span. For the exception of the selected experimental parameters, all sources of stochasticity (e.g., assignment of agents' price prediction models, learning time horizon, income) are held constant using a random number seed. .
The model is initialized with a developed CBD that contains a mix of housing types on lots of either . , . , , or acres and house sizes of , or , square feet ( Figure ). All new houses that are built over the course of model execution are a combination of one of these existing house and lot sizes (i.e., no new housing types are introduced). At the beginning of the simulation, the same number of consumers are initialized as there are number of houses in the CBD. Consumers are heterogeneous in their incomes and preferences for house size, lot size, and proximity to the coast. Consumer incomes and housing and location preferences are randomly drawn from uniform distributions. Undeveloped land is subdivided into -acre parcels ( by cells) that are regularly distributed across the landscape (Figure ). Each undeveloped parcel is assigned to a single landowner agent. ). The developer forms a bid price for each undeveloped parcel based on spatially explicit expected rents for housing, and landowners form asking prices based on their expected value of the land in its most profitable use (e.g., agriculture versus residential housing). If the developer's bid price is higher than the landowner's asking price, then a transaction occurs and the entire parcel transfers ownership. .
Given observed land purchase prices, the developer ranks each housing type from highest to lowest expected return (i.e., expected rents net of construction and land costs) given distance to the CBD and coast in every owned and/or undeveloped cell in existing land holdings or newly acquired parcels. The housing type with the highest expected return in a given cell is built first, with others following in remaining owned and undeveloped cells in the order of descending expected returns. Housing construction of a given type continues until as many houses of that type are built as the developer expects will be demanded by consumers, or all vacant land is developed. New houses are assigned the closet distance from the CBD and highest coastal amenity value of any cell in the new lot's footprint. New and existing vacant houses are then o ered for sale in the housing market with an asking price equal to the developer's expected rent for the given housing type and location.
. New housing consumers are introduced into the 'consumer pool' of existing consumers at the start of each time step at a rate of ten percent a year. Each consumer c generates utility from housing and a composite nonhousing good, x. Each house i is characterized by its size, h i , its lot size, l i , and an amenity that is a function of the distance from the coast, a(d i ), where d i is house i's distance from the coast. We assume that the utility function has a standard Cobb-Douglas structure: . A consumer has income, I, pays (annualized) price, P i , for house i, and pays travel costs, ψ i , which depend on house i's distance from the central business district (CBD). Consumers also face an annual probability, ρ, of property damage from coastal storms. When a storm occurs, it imposes a cost, which we assume is a percentage of the house's value. The percentage varies with the house's distance from the coast, C(d i ); where δCi δdi < 0 Thus, the budget constraint is given by: In this version of the model, we assume that consumers know the probability of a storm and the costs they will incur if a storm takes place. We also abstract from considerations of insurance. Both assumptions will be relaxed as a focus of future research. .
Consumers place bids on all houses that generate positive utility and are thus within their budget constraints given asking prices. Each consumer observes the number of other consumers placing bids on their set of prospective houses and adjusts their bid according to competition levels. If there are more bids than there are houses being bid on, then bid prices are adjusted upwards. If there are fewer bids than there are houses being bid on, then bid prices are adjusted downwards. Bid adjustments are designed to reflect the relative supply of and demand for housing and quantify competition with the housing market. Houses are assigned to consumers with the highest bids for each house. If a consumer has the winning bid on multiple houses, then the consumer chooses the house that generates the highest utility given the winning bid prices. Consumers that locate in a house are assigned a residence time drawn randomly from a normal distribution (µ = 12.5, σ = 11 time steps). When residence time is exceeded, the consumer moves-out, re-enters the consumer pool, and the newly vacant house is put back on the market. Consumers that do not locate a er three consecutive time steps are removed from the consumer pool.

Space-and Time-Varying Global Sensitivity Analysis
. The overall objective of the following GSA model experiments is to elucidate the extent to which parameters and their interactions contribute to variability in model outcomes. These experiments do not to establish causality between any particular parameter and development patterns in any particular location or single model realization. The intent is to illustrate the utility of combining variance-and density-based GSA approaches, rather than to draw conclusions about any particular coastal systems.
Variance-based GSA decomposes the variance (V ) of a selected model output (Y ) from the perturbations of k = 7 model inputs ( Table ) into the contributions of single (V i ) and combinations of model inputs with increasing dimensionality (Ligmann-Zielinska & Sun ): where V i is the contribution of model input X i to the overall variability in model output, and V ij , for example, is the variance in model output due to interactions between model inputs X i and X j . Variance decomposed thusly is then used to calculate first-order (S i ) and total-e ect (ST i ) indices for every model input i in i = 1 : k: The proportion of model output variance (V ) attributed to input i independent of other inputs (k − 1) is given by S i , and ST i describes the overall contribution of input i and its interactions with other inputs to the overall model output variance (Ligmann-Zielinska & Sun ; Saltelli et al. ). The contributions to variance in model output Y of all inputs other than i is represented by the expression

Description PDF Statistics
Storm Frequency (M f ) Probability of a coastal storm of a fixed severity occurring in any given year. Economic discount rate used to annualize land and housing prices.
N ( . , . ) Table : Description of selected uncertain model input parameters and their distributions used for quasirandom sampling (refer to text for explanation). TEYhe probability density function (PDF) column specifies whether the distribution used was discrete (D), continuous and normal (N ), or continuous and uniform (U ). Descriptive statistics specify the minimum and maximum values (D or U ), or mean and standard deviation (N ).
. We follow the time-varying GSA method developed by Ligmann-Zielinska & Sun ( ), which is based on the quasi-random Sobol estimation procedure (Saltelli ; Lilburne & Tarantola ) and implemented using a sequential estimation algorithm in Matlab developed by Tarantola ( ). In addition, building on the approaches by Hall et al. ( ), Lilburne & Tarantola ( ), and Ligmann-Zielinska et al. ( a), spatially distributed sources of sensitivity are investigated at the regional (aggregate) scale as well as in several representative sub-regions within the landscape (fine-scale). For time-GSA, S i and ST i are calculated for every time step to determine to what extent particular parameters are stronger or weaker drivers of model outcomes over time. For space-GSA, the landscape is divided into five 'zones' (i.e., sub-regions) based on their perceived value for development relative to their proximity to the coast (Figure ). In preliminary experimentation, each zone showed internally consistent yet distinct development dynamics from other landscape regions. The CBD is excluded from analysis because the selected model outputs -number of houses built, percent developed area, and diversity of the housing stock -remain constant over the course of model execution in that location.
The density-based GSA approach developed by Pianosi & Wagener ( ), PAWN, is both methodologically and conceptually complementary to the variance-based GSA methods. PAWN is based on the comparison of conditional and unconditional CDFs of model output. The unconditional CDF is simply the distribution of model outputs generated from all parameter values used in the Sobol estimation procedure for the variance-based GSA. The conditional CDF is generated by fixing the value of one input parameter while letting other parameter values vary. Fortunately, this approach can be applied to existing datasets (i.e., using the same simulation data for both variance-and density-based analyses) by grouping model executions with similar input parameter values, which can then be treated as 'fixed' value sub-sets from the parameter distribution (Pianosi & Wagener ). Conditional distributions account for changes in outcome distribution when variability from a particular parameter (or parameter group) is removed, and thus the distance between the conditional and unconditional distributions is a measure of output sensitivity. The significance of the distance between conditional and unconditional distributions is tested with the Kolmogorov-Smirnov (KS) statistic (Kolmogorov ; Smirnov ). .
Fixed groups are specified here by sub-setting each input parameter distribution generated by the quasi-random Sobol estimation procedure into ten equal-frequency quantiles, and thus the 'fixed' value for each quantile group is the median of the quantile. Conditional CDFs are then created by removing one parameter group at a time to observe the e ect on model outcome distributions. All conditional CDFs are then plotted against the unconditional CDF and the KS statistic is calculated for each fixed value and compared against the critical value of the KS statistic at confidence level of . . KS statistic values exceeding the critical value indicate a statistically significant influence of the parameter at the given fixed value. Importantly, this analysis is applied only to final model outcomes, as no spatial or temporal details about parameter sensitivities are provided. Thus, PAWN is a complement to variance-based in that it provides insight into where in parameter distributions sensitivities in model outcomes are produced.

Experimental approach .
A key consideration when conducting GSA is the number of model executions required to obtain accurate sensitivity estimations. Estimations of sensitivity indices S i and ST i are calculated from a number of base samples (N ), which each contain k + 2 quasi-random parameter samples used for sequential estimation, and bring the total number of model executions required to N (k + 2) (Tarantola ). The larger the N, the more precise the sensitivity estimates, but practically the selection of N depends on the computational cost of the model. Relatively inexpensive models can have an N in excess of , but more expensive models may have an N on the order of - (Tang et al. ). Additionally, as the complexity and nonlinearity of a given model increases, so too should N . For this study, an N = 100 is selected for a total of model executions as a compromise between estimation accuracy and computational demand. A single model execution using the parallel computing toolbox in Matlab requires roughly minutes run-time, which results in a total of , hours of run-time on a multi-threaded . GHz AMD Opteron™ Processor .

Input parameter distributions .
Seven independent model input parameters are used (Table ): storm frequency, storm severity, maximum natural amenity value at the coast (A o ; Figure ), rate of coastal amenity decrease with distance from the coast (r; Figure ), range of consumer preference for coastal amenity (δ c , Equation ), travel costs to CBD (ψ i , Equation ) , and economic discount rate. The economic discount rate is used by the developer to calculate the profitability and formulate future asking prices of particular housing types given observed housing and land prices. . For all normally distributed inputs we assume a standard deviation of half the mean to ensure su icient perturbations of parameter values. Probabilities of storm occurrence (M f ) are aggregated from Costanza et al. ( ) across all storm severity categories for the Atlantic coast to estimate discrete annual probabilities of storm occurrence of any severity. Storm severity (M s ) is varied independently by modifying percent loss of property value (Figure ). Storm frequency and severity are decoupled in this way for use with an expected cost framework for resolving consumer location decisions. Consumer preferences for coastal amenities (δ c , Equation ) are randomly drawn from a uniform distribution of varying range. For the purposes of exploring model output sensitivities to this formalization, the quasi-random sampling procedure is applied to the maximum of the preference distribution (e.g. [ . , . ] vs. [ . , . ]) to explore the e ects of a more or less heterogeneous and strong preference for coastal amenities among the consumer agent population.

Output parameter selection .
Three model outputs are selected for sensitivity analyses: number of houses, percent developed area, and diversity of the housing stock. The first two outputs relate mainly to spatial patterns of development and provide measures of the cumulative e ects of household location decisions over time. The diversity of the housing stock, or housing stock evenness, at each period is measured using Shannon's entropy. This index describes the relative likelihood of a particular house type being built in a given period. In other words, as the distribution of housing types in the overall housing stock becomes more even (Shannon entropy approaches ), each housing type is equally likely as any other to be built and housing o erings are essentially random. Conversely, as the distribution of housing types in the overall housing stock becomes more uneven (Shannon entropy approached ), a particular housing type (or set of housing types) is more likely to be o ered reflecting relatively more desirable and/or profitable housing types. Such a measure can be related to self-reinforcing and/or pathdependent processes in the development process.

Results
. The complete set of model outputs included temporally and spatially explicit land-use maps giving the locations, types, and prices of housing, consumer preferences and incomes, and transaction prices for land. Selected model outputs were aggregated within each of the sub-regional 'zones' (Figure ) and sensitivity indices were calculated every time step for each zone. Figure shows an example of variation in development patterns with low and high coastal amenity decay rates (r) for two selected model executions. Note that a perturbation in only this input parameter resulted in large variation in the location of development.

Sequential sensitivity estimation errors .
The accuracy of sensitivity indices depends on the selection of N (Saltelli et al. ), however the optimum N is o en not known a priori for complex models. Using the quasi-random Sobol sampling method and algorithm for sequential estimation of sensitivity indices developed by Tarantola ( ), the root mean squared error (RMSE) of first-order sensitivity was calculated for each model execution to track convergence of sensitivity estimates to a moving average for each selected model output (Figure ). Given the computation expense of the coastal CHALMS model, RMSE values below percent (RMSE = . ) were deemed satisfactory (Saltelli et al. ). RMSE values for the number of houses built and percent developed area (Figure a and b) for all input parameters were generally less than the threshold for all model executions, except for slight increases in variability due to discount rate. In the case of housing stock evenness (Figure c), variability due to storm frequency and discount rate were above the RMSE threshold by the final model execution. However, RMSE values were relatively stable, fluctuating around their overall mean values through all model executions. Hence, N = 100 samples were su icient to conduct both GSAs on these model outputs. For variance-based GSA, time-dependent first-order (S i ) and total e ect (ST i ) sensitivity indices are plotted cumulatively at the regional and zonal levels for each selected model output. These are paired with outcome variance from density-based GSA presented as CDF and KS statistic plots of regional outputs from input parameters. For interpretation of variance-based sensitivity indicators, the e ects of model inputs are additive when interactions are linear and the sum of S i nears (Ligmann-Zielinska & Sun ). As non-additive behavior increases, the proportion of variance attributed to interactions among parameters increases, the sum of Si decreases, and areas plotted as white in Si plots increase. Regionally, variance in houses built and percent developed area was mostly additive, while variance in housing stock evenness was mostly due to parameter interactions. For all outputs, no single input parameter was consistently influential over time, as sensitivity indices oscillated frequently as the development process unfolded. Furthermore, sensitivity profiles at the zonal level were quite di erent from one another and regional profiles, which emphasized the importance of studying spatial as well as temporal global sensitivity.

Number of houses built .
At the regional level, storm frequency (M f ) was responsible for much of the variance (S i ) observed in the number of houses built (Figure ) a er seventh model time step. On average, amenity preference (δ c ), travel costs (ψ i ), and discount rate (R) accounted for an average of % of total variance (ST i ) when interacting with other factors. Closer inspection with density-based GSA shows that outcome variance from storm frequency at the regional level is driven exclusively by the highest probability storm climate (Figure a). Additionally, the discount rate consistently influenced the number of houses built, with high and low discount rates and moderate discount rates tending to lead to fewer and more houses than the global average, respectively (Figure b). . Zonal sensitivities were quite di erent than those at the regional level. Discount rate (R) and travel costs (ψ i ) were almost entirely responsible for variance in the number of houses built in zone . Storm frequency was not influential in any zone, whereas the other input factors alternated in importance over time. Inputs amenity preference, amenity decay rate, travel costs, and discount rate are all influential factors in consumers' decisions to locate in proximity to either the CBD or coast, as well as interacting strongly to drive housing and land prices. For example, competition for coastal property in zone may have pushed consumers inland and driven development in zone .

Percent developed area
. At the regional level, storm frequency was responsible for much of the variance observed in the percent of developed area (Figure ) a er seventh model time step. Just as for the number of houses built, density-based GSA revealed that the highest storm frequency was mostly responsible for regional variance. Less area developed on average with high and low discount rates, while moderate discount rates lead to more developed area than the global average. Figure a and b also show that higher values of amenity preference and decay rate influenced the percent developed area, with increased preferences for amenities and faster decay of amenity values led to less developed area overall. However, because the relative di erences between STi and Si were large at various time steps for model inputs other than storm frequency, we conclude that all contributed significantly to overall variance through their interactions with other model parameters (i.e., non-additive; Ligmann-Zielinska & Sun ).

Figure :
Time-varying first-order (S i ) and total e ect (ST i ) sensitivity indices for the percent developed area are color-coded for the contribution of each model input. Regional results are plotted in the first row, while zonal results are plotted on the lower two rows.
. Zonal sensitivities were quite di erent from those at the regional level, as no single model input was consistently important over time across all zones. All zones exhibited cyclical patterns of influence on variance with one or two parameters dominating in any single time step. In general, travel costs and discount rate consistently accounted for % to % of variance, while amenity-and storm-related decreased and increased in importance, respectively, moving from inland to coastal zones. Unlike the number of houses, storm frequency remained important at the zonal level for the percent developed area. This suggests that increased storm frequency tends to modify the density of housing, because variance in developed area is high and the number of houses is comparatively insensitive to storm frequency changes.

Housing stock evenness .
Sensitivity of housing stock evenness was mainly driven by interactions among inputs at both regional and zonal levels (Figure ). The sum of Si values averaged % throughout model execution at the regional level. Similarly, ST i values were relatively evenly distributed among all of the model inputs further reinforcing the importance of input interactions for housing stock outcomes. Density-based GSA was particularly revealing in this context. Higher amenity preference and lower amenity decay rate values tended to create more even housing stocks (Figure a and b). Both of these conditions provide more overall amenity value across the landscape, which mediated the e ects of income and location di erences that tend to favor some housing types over others. For example, ceterus paribus, higher income consumers tend to preferentially occupy houses on larger lots, and high density housing tends to occur in close proximity to the coast due to higher profit margins. Discount rates also demonstrated directional influences of housing stock evenness, with lower rates leading to more even housing stock (Figure c).

Discussion
. Space-and time-varying GSA revealed input parameter interactions di ered over both time and spatial scales and across selected model outputs. For the exception of housing stock evenness, the sum of first-order e ects at the regional scale ranged from % to %, meaning that a high proportion of output variance at that scale was additive. However, there were instances in particular zones and at various points during model execution that the sum of first-order e ects dropped to -% across all model outputs, indicating the importance of localized spatial interactions in the development process. For example, storm frequency was highly interactive with other inputs in non-coastal zones, specifically in Figures and , which could be indicative of spatial 'spillover' e ects from high land and/or housing prices in the coastal area (Irwin & Bockstael ). In other words, market dynamics in one part of the landscape influenced development in adjacent areas of the landscape, i.e., 'spill-over'. .
Furthermore, input interactions over time were apparent in the cyclical patterns of variance observed at both regional and zonal levels. Cyclical patterns were observed in total-e ect indices for all regional level outputs as inputs alternated in importance of driving output variance over time. Two intrinsic features of the model that reflect land and housing markets characteristics explain the emergence of these patterns. Land markets have been described as "thin markets" (Irwin & Wrenn ), meaning market transactions occur at relatively low frequency and are geographically dispersed, as opposed to other markets (e.g., financial) which tend to have many more transactions and more continuous price dynamics as a result. Additionally, the discrete nature of land parcels and housing options over which agents make decisions create discontinuous land and housing market dynamics. Such cyclical behavior is o en attributed to delayed feedbacks created by a time lag between a change and the e ects of that change on the rest of the system (Manson ). This form of complexity is important to represent in ABMs, particularly when investigating spatially-explicit interactions that operate unevenly over time across the landscape. The relative regularity of the cyclic patterns is explained by the constant, imposed rate of population growth and fixed landowner parcel sizes. In other words, demand for land and amount of land purchased to meet that demand were relatively steady over time, which produced a consistent (although not deterministic) cyclic structure in model sensitivity. Time-varying GSA enabled the detection of these interactions via their e ects on output variation.

.
Sensitivities of housing stock evenness, in particular, illustrated the highly complex and nonlinear nature of market processes. Housing stock diversity results from many short timescale feedbacks between structural features of land and housing markets (e.g., transaction costs, construction costs) and behavioral features of landowners, developers, and housing consumers (e.g., price expectations, location preferences, risk tolerances). Sensitivity indices showed that the majority of variance in housing stock diversity was due to interactions among input parameters and model features (Figure ). Application of conventional OAT approaches to sensitivity analysis would not have detected the importance of input interactions in such emergent phenomenon, and may have even generated misleading conclusions about the magnitude and relative importance of housing stock evenness sensitivity to particular input parameters. .
The relatively large di erences in sensitivities across zones, coupled with evidence of the importance of interactions for explaining variance, suggest that housing stock diversity is potentially a locally path-dependent and regionally state-dependent phenomenon. Page ( , p. ) defines a state-dependent process as one in which the outcome at any time depends only of the state of that processes at a given time, and di erent states that lead to the same future state are equivalent. In this context, the relative proportion of di erent housing types present in the overall housing stock at any given time, described here as the evenness of the housing stock using Shannon entropy, can be considered the housing stock state. A path-dependent process is one in which the outcome at any time depends on both past states and their order of occurrence Page ( , p. ). At the zonal level, housing stock evenness was generally the result of path-dependent processes. Local spatial interactions and market forces made it more likely that a particular housing type would be built in each zone if it had been built there previously (Figure ). In contrast, regional level housing stock evenness tended to become more even over time, suggesting that regional housing stock diversity was the result of state-dependent processes. Feedbacks from the regional land and housing markets stabilized overall housing diversity over time and ensured that roughly similar housing stocks consistently emerged at the regional level regardless of the order in which particular housing types were built in di erent locations. Such complex, cross-scale dynamics would not have been detectable with conventional sensitivity analyses. .
While space-and time-varying GSA is an appropriate analytical tool for complex ABMs, its application presents a number of conceptual and computational challenges that conventional methods do not. The number of model executions required for GSA can be exponentially higher than for convention methods. Su icient sample size is critical to producing accurate estimates, with more uncertain outputs potentially requiring (many) more than other outputs. Given the computational expense, the analyst is faced with a trade-o between extended processing time and estimation accuracy. Sequential estimation techniques, such as the Sobol method used here (Tarantola ), have the potential to optimize the number of model executions performed (Tang et al.  ). However, this method prohibits parallel model executions, which can lead to exorbitant run times for computationally expensive models. Here, we selected a large number of initial base samples based on the literature (Tang et al. ; Ligmann-Zielinska & Sun ; Tarantola ), generated quasi-random parameter samples beforehand, executed each model run in parallel, and evaluated convergence of sensitivity estimates post-hoc (Section . ). This approach required , hours of run-time as opposed to approximately , hours that would have been required if model execution had been performed in sequence.

.
One of the benefits of GSA is the identification of uncertain input parameters that contribute insignificantly to output variance, and thus can be fixed without changing overall model output (Ligmann-Zielinska & Sun ; Ligmann-Zielinska et al. a). However, as this analysis demonstrated, in highly nonlinear ABMs, it may not be possible to identify and fix such parameters because of their highly interactive nature. If particular parameters can be specified empirically, such as travel costs and discount rate, this may eliminate sources of uncertainty even though the model is sensitive to these parameters. In contrast, some input parameters to which the model is highly sensitive, such as the spatial distribution of amenity values and consumers' location preferences for coastal amenities, are likely to be heterogeneous across consumers and di icult to estimate empirically. In such cases, GSA is a useful tool for describing model sensitivities even if it does not ultimately lead to reducing model sensitivities through parameter fixing. In either case, GSA can identify parameters responsible for significant output variance, which then become priorities for data collection when developing an empirically parameterized model.

Conclusions
. Although the CHALMS model is initialized with stylized parameter values and landscape configuration, this application of combined variance-and density-based GSA highlighted issues that apply to both theoretically-and empirically-based spatial ABMs. Variance-based GSA can be extremely computationally expensive, particularly when model outcomes are highly uncertain and require thousands of model executions to obtain accurate sensitivity estimates. As we demonstrated here, an iterative process in which quasi-random input parameters are generated a priori, model executions are performed in parallel, and estimation error is calculated post-hoc is a promising alternative to sequential model execution and sensitivity estimation for models with equally or more complex behavior than CHALMS. Density-based GSA using PAWN (Pianosi & Wagener ) also provides insight into the magnitude and directionality of model sensitivities, and it can be applied to existing model data, eliminating the need for additional model executions. .
The computational cost of GSA methods for complex models is daunting but necessary. Once a model is capable of producing nonlinear and adaptive behaviors, such as path-dependent development patterns, a threshold is passed beyond which model outcome variability, and therefore the number of model executions required for GSA, is not reducible (Brown et al. ; Manson ). Such a gradual approach is both pragmatic, given the computational demands required, as well as necessary in order to implement comprehensive sensitivity analysis of complex ABMs.

.
Spatially and temporally explicit variance-based GSA and complementary density-based GSA have the potential to deepen understanding of the sensitivities of complex, dynamic ABMs. Uncertainty in input parameters contributes not only to sensitivity of model outcomes, but also to spatially-and temporally-varying sensitivities in model dynamics. One of the main advantages of ABMs over other top-down modeling or statistical approaches, particularly for studying land-use change, is to explicitly model the time path, or evolution, of emergent phenomenon. However, interpretation of results remains more challenging than for conventional methods, especially when ABMs possess many nonlinear characteristics, such as path-dependence (Brown et al. ; Manson ). This emphasizes the need for parsimony in model design. Even simple models can produce complex behavior (Sun et al. ), but the task of GSA is made easier and more e icient if the number of uncertain parameters (or overall number of parameters) can be reduced. .
The ABM research community is engaged in an ongoing debate about the most appropriate SA methods. Specifically, Saltelli et al. ( ) and Ligmann-Zielinska et al. ( a) have concluded that variance-based GSA is the most appropriate method for process-based models characterized by parameter interactions and emergence. Conversely, Ten Broeke et al. ( ), recommend OAT over regression-based and GSA methods, because, they argue, these other methods consider only averaged parameter e ects, masking nonlinear model mechanics.
In light of our findings here, we agree with ten Broeke et al. on this point. Variance-and density-based GSA methods used alone are not su icient for understanding model mechanics, because they fail to address either the directionality or dynamics of parameter e ects on model outcomes, respectively. However, also clear from our findings is that accounting for parameter interactions is essential for understanding and quantifying the full scope of model behavior driven by uncertain parameters. Because OAT methods do not consider interaction e ects, they will always be limited to partial examinations of model outcome space. Emergent model behaviors resulting from unique, simultaneous combinations of two or more parameters cannot be captured by perturbing one parameter at a time. Only global approaches that consider the full range of all uncertain parameters within a single analysis can rigorously investigate such behaviors. Furthermore, we have demonstrated that pairing complementary GSA methods can address the remaining criticisms of GSA methods. Specifically, density-based GSA can reveal which parts of individual parameter distributions are associated with model output sensitivities, and how robust model outcomes are to various parameter combinations within and across multiple parameter distributions (i.e., hotspots within the global parameter space). The prevalence (or sometimes dominance) of parameter interactions in ABMs requires the comprehensiveness that GSA methods provide, and restricts OAT methods to a useful means of model exploration. .
The unique combination of GSA methods implemented here enriched our understanding of how our assumptions about parameters values a ected modeled dynamics. Highly sensitive model outputs were identified for future model experiments in which the sampling space will be refined, and data collection e orts will be targeted to constrain parameters that produce model sensitivity. Critical sources of model sensitivities were also identified across space and over time that lead to new research questions. For example, the influence of storm frequency was stronger in certain locations and at di erent times during simulation runs. This leads to the question: how does the timing of storms -earlier or later in simulations, which impact less or more developed area, respectively -a ect development dynamics? Ultimately, these are the types of questions that comprehensive GSA can raise, and force us as modelers to deepen our understanding of complexity in land-use ABMs.