Two-Step Sensitivity Testing of Parametrized and Regionalized Life

May 8, 2013 - Land Use in Life Cycle Assessment: Global Characterization Factors Based on Regional and Global Potential Species Extinction ... Ignorin...
0 downloads 0 Views 859KB Size
Article pubs.acs.org/est

Two-Step Sensitivity Testing of Parametrized and Regionalized Life Cycle Assessments: Methodology and Case Study Christopher L. Mutel,†,* Laura de Baan,‡ and Stefanie Hellweg† †

ETH Zurich, Institute of Environmental Engineering, 8093 Zurich, Switzerland ETH Zurich, Institute for Environmental Decisions, 8092 Zurich, Switzerland



S Supporting Information *

ABSTRACT: Comprehensive sensitivity analysis is a significant tool to interpret and improve life cycle assessment (LCA) models, but is rarely performed. Sensitivity analysis will increase in importance as inventory databases become regionalized, increasing the number of system parameters, and parametrized, adding complexity through variables and nonlinear formulas. We propose and implement a new two-step approach to sensitivity analysis. First, we identify parameters with high global sensitivities for further examination and analysis with a screening step, the method of elementary effects. Second, the more computationally intensive contribution to variance test is used to quantify the relative importance of these parameters. The twostep sensitivity test is illustrated on a regionalized, nonlinear case study of the biodiversity impacts from land use of cocoa production, including a worldwide cocoa products trade model. Our simplified trade model can be used for transformable commodities where one is assessing market shares that vary over time. In the case study, the highly uncertain characterization factors for the Ivory Coast and Ghana contributed more than 50% of variance for almost all countries and years examined. The two-step sensitivity test allows for the interpretation, understanding, and improvement of large, complex, and nonlinear LCA systems.



INTRODUCTION Sensitivity analysis is part of the iterative process of improving life cycle assessment (LCA) model inputs to decrease result uncertainty and provide decision support. Global sensitivity analysis, however, is rarely found in published LCA studies, and an approach that is comprehensive, understandable, and integrated in LCA software would facilitate routine sensitivity analysis. The term “sensitivity analysis” is used in different ways in the LCA literature, including the identification of product system improvement potentials, but in this manuscript we follow the definition of Saltelli et al.: “The study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input.”1 This differs from uncertainty analysis, which quantifies the uncertainty of model outputs. Four main types of self-identified sensitivity analysis can be found in the LCA literature. The most popular is one at a time (OAT) variation, also called perturbation analysis.2,3 In OAT testing, the practitioner manually varies one parameter by a certain percentage and computes the change in model outcome, calculating the sensitivity factor of that parameter. A variant of OAT testing is scenario analysis, where scenarios have different values for several parameters. Although OAT approaches meet the ISO LCA standards,4 OAT testing has been called “perfunctory sensitivity analysis”, as it ignores interaction effects between parameters, does not assess all model © XXXX American Chemical Society

parameters, and assumes all parameters have strictly linear effects.5 OAT testing also usually varies each parameter by a set percentage, ignoring the probability distribution and expected ranges of each parameter. Both manual OAT testing and scenario analysis furthermore rely on the practitioner to select the relevant parameters for variation. However, the relevant parameters are often not known, and sensitivity analysis is specifically used to identify the most important varying parameters.6 Second, the structure of a matrix-based LCA can be used to calculate sensitivity factors for all parameters.7,8 This second type of sensitivity analysis examines all system parameters, but does not include the expected range of each parameter, and hence ignores the actual contribution of each parameter to model variance. However, the formulas needed are relatively complex and assume a certain number and order of matrices. For application to regionalized LCA, the formulas would need to be updated,9 and they cannot be applied to parametrized databases, which include nonmatrix calculations. A third approach is to use the variance in model parameters in the matrix-based calculation of sensitivity factors.10 This Received: December 12, 2012 Revised: May 8, 2013 Accepted: May 8, 2013

A

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

calculated from a set of Monte Carlo iterations. The Spearman’s rank-order correlation coefficient (ROCC) of each parameter’s values with the LCIA scores for the set of Monte Carlo iterations is calculated. Rank-order correlations can range from −1 to 1. These ROCCs are then squared, which ensures that the resulting values are positive, and gives proportionally more weight to higher correlations. The final CTV values are each parameter’s squared ROCC divided by the sum of the squared ROCCs for all model parameters,11 as shown in eq 1:

approach calculates the propagation of error using the variance of each parameter probability distribution. This sensitivity analysis can be quickly computed, but suffers from the same drawbacks as the earlier matrix-based techniques. In addition, using only the variance of a probability distribution can neglect significant information for asymmetric or empirical distributions. Finally, contribution to variance (CTV) has been proposed as a global sensitivity test for life cycle assessment.11,12 CTV, as the name implies, expresses the contribution of each parameter to the overall variance in model outcomes. This relative contribution is especially helpful in deciding where to improve data quality to decrease overall model variability. In contrast with analytical methods, CTV is nonparametricit does not assume any particular parameter distributions or mathematical structure of the model, but operates only on the Monte Carlo sampled model inputs and calculated outputs. It is also easy to calculate and interpret. CTV has been applied to a small set of input variables, but has not yet been successfully applied to entire inventory databases. The next generation of LCA databases, such as ecoinvent version 3, will include both regionalization and parametrization,13,14 and robust sensitivity testing will be crucial for understanding the resulting large and complex systems. Regionalization will dramatically increase the number of parameters in the model, and will increase the importance of considering impact assessment uncertainty, as characterization factors will differ from region to region. For example, one regionalized method for assessing water use provides characterization factors for over 10 000 watersheds,15 and we have modeled over 5000 electrical generating units in the United States alone.9 Parameterizationthe ability to use variables and equations in life cycle inventoriescan introduce nonlinearities. In our case study of globally traded cocoa products, crop yield is expressed as production divided by cultivated area, and therefore area is not linearly related to yield. Sensitivity tests need to consider not just the traditional parameters in the biosphere, technosphere, and characterization matrices, but also these new parametrized inventory variables. In this manuscript, we address several impediments to the routine use of comprehensive sensitivity analysis in LCA. To understand the computational limits of CTV, which is one reason why CTV has not been applied to entire inventory databases, we use numerical simulation to determine the number of Monte Carlo iterations needed for CTV values to converge. We then propose a two-step sensitivity analysis that can be applied to very large LCA models, and which has several advantages. To facilitate the inclusion of this two-step approach in LCA software, we provide an open source implementation of our method. We illustrate this new sensitivity analysis on a case study of the biodiversity impacts of land use occupation impacts from cocoa production with over 30 000 uncertain parameters. The case study includes a large number of nonlinear parameters, including regional characterization factors modeled with kernel density estimators. The case study model also proposes a simplified approach for modeling transformable trade commodities in the absence of specific data on domestic value chains, stocks, or consumption.

CTVi =

ROCCi 2 ∑j ∈ n ROCCj 2

(1)

where i is the parameter whose CTV is being calculated, and n is the set of all stochastic parameters (including i) in the model. The formal mathematical properties of Spearman’s rankorder correlation have been thoroughly investigated.16,17 However, no simple equation exists for what we will term the expected uncorrelated error: the error due to nonzero ROCCs for parameters which have no correlation with the model output. To determine the expected uncorrelated error, we numerically simulated a set of random variables, only one of which was set to be the model result, and therefore was perfectly correlated. The Supporting Information (SI) gives details on how, by varying the sample size, type, and number of variables, we could investigate the mathematical properties of CTV. Method of Elementary Effects. Because CTV is relatively computationally inefficient, we would like to have a global sensitivity method that could be used to screen for sensitive parameters. The method of elementary effects (MoEE) is a robust method for the identification of global sensitivity and nonlinearity of model parameters,18,19 but has not yet been widely applied in LCA studies. De Koning et al.20 focused on the Fourier amplitude sensitivity test, but also applied the MoEE in their SI. A detailed explanation of the MoEE is provided in Saltelli et al.,21 but we examine the basic concepts. The MoEE examines the sensitivity of model outputs to changes in each parameter, but instead of keeping the other parameters at their median or mean values, it systematically varies all parameters in series to explore the full range of model outcomes. Uncertain parameters should not be varied by a uniform percentage, as this ignores the different expected ranges of each parameter. For example, the amount of fossil fuel combusted in an electrical generator per unit of energy has a relatively narrow range, while the toxic impacts from the combustion emissions could vary by orders of magnitude. We want to instead compare the same percentile change in each parameter’s range. To use the MoEE, therefore, all uncertain parameter values are first linearly scaled from their actual uncertainty ranges to the (0, 1) unit interval. In most cases, this can be accomplished by a subtraction of the minimum value and a subsequent division by the maximum value. The SI shows how the cumulative distribution function can be used to transform data values to percentages. This translation allows for the same relative change for each parameter in their expected ranges, and therefore for the direct comparison of parameter sensitivities. A large number of trajectories are constructed, and each point in a trajectory represents a set of parameter values. In each trajectory step, only one parameter value varies. Each



MATERIALS AND METHODS Contribution to Variance. We begin by trying to understand the computational reasons why CTV has not been successfully applied to large inventory databases. CTV is B

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

practice, it is easy to include more parameters in secondary sensitivity analyses, if needed. Case Study. Case Study Overview. We applied the twostep sensitivity testing approach to a case study of the land use occupation impacts on biodiversity of cocoa production. Cocoa was chosen for the case study because it has a relatively short value chain, a limited number of useful byproducts, large production and trade volumes, and is produced in areas with high potential biodiversity. We modeled the global production and trade of cocoa and cocoa products from 1997 to 2007, and assessed one kilogram of cocoa-equivalent from country market mixes. The coupled production, trade, and value chain model includes uncertainty information for all parameters. Different techniques to model the uncertainty distributions were used for each parameter type based on the available data. For the impact assessment, we developed regional cocoa-specific land use characterization factors and their respective uncertainty distributions. In addition to understanding the differences in the impact of cocoa consumption in various countries, we used the gathered data to assess changes in trade patterns and impact through time. Case Study Model Construction. Cocoa is contained in different concentrations in a range of traded commodities. As we did not have commodity-specific data on domestic consumption or stocks, or on the domestic value chain of different cocoa commodities, we constructed a simplified market model, trading “cocoa equivalents” according to the average cocoa mass in each commodity. The computational model was similar to standard matrix-based LCA calculations, with a vector of characterization factors, a vector of biosphere flows, a demand vector, and a technosphere matrix. Biosphere flows are a vector because only one environmental intervention, land occupation, is calculated. The technosphere matrix has an upper and lower half, with rows in the top half representing the fractional share of domestic production of cocoa in each country, if any, in that country’s cocoa market mix, and rows in the bottom half representing the fractional share of imports from the corresponding export country into the importing country’s market mix. The sum of all imports and trade production in each country’s market mix is one. The technosphere matrix has twice as many rows and columns as the total number of countries considered; in the case study, the technosphere matrix had 534 rows and columns. The SI has graphics showing the construction of the characterization and biosphere vectors, and the technosphere matrix. The case study mathematical framework was

trajectory should have n + 1 steps, where n is the total number of parameters. SI Figure S1 shows an example of six trajectories of three parameters each. Although SI Figure S1 has only three parameters, one can imagine a hypercube of many more parameters, with the chosen trajectories acting similar to spacefilling curves. In each trajectory step, one parameter value will change between four preselected values, from the zeroth to the 66th percentile of its distribution, or from the 33rd to the 100th percentile, or vice versa.19 In the case study, we combined uncertainty distributions with defined limits, such as the uniform distribution, and uncertainty distributions without such limits, such as the normal distribution. There is no sensible value for the zeroth percentile of a normally distributed parameter, so it was necessary to choose a restricted percentile range. To be consistent, we therefore used the 2.5, 34, 66, and 97.5 percentiles for all distributions. The discerning power of MoEE does not appear to increase when more percentile steps are used.19 From this set of candidate trajectories, a group of trajectories are chosen. The Kennard−Stone algorithm is used to pick the group of trajectories that approximately maximizes the Manhattan distancethe sum of the absolute differences over each axisfrom each trajectory to all other trajectories. Common practice is to choose 10 trajectories from 1000 candidates.19 The model is evaluated for each point in each trajectory. Parameter values should be transformed back from the (0,1) unit interval to the true parameter values before model evaluation. There are 10 independent differences in model outcomes due a change in each parameter values, and one can calculate the average and standard deviation of these 10 responses. The μ* statistic is calculated from the absolute value of the model responses to each parameter. μ* is a more robust indicator than the average because it assesses only the magnitude of the change, not the sign, and is not susceptible to averaging effects when the response can be both positive and negative.19 The μ* statistic is a measure of global sensitivity, and the standard deviation, σ, is a measure of nonlinearity of parameter changes relative to model outputs, as it measures the differences in model output to a fixed relative change in model inputs. The MoEE computer code we developed is open source and is described in the SI. A New Two-Step Approach for Sensitivity Testing. We propose a two-step method for sensitivity testing of LCA system models that combines CTV and the MoEE. First, the MoEE is used to calculate the global and nonlinear sensitivities of all parameters in the model, including impact assessment characterization factors. Parameter sensitivities are ranked, and parameters with high global or nonlinear sensitivities are chosen for further analysis. Contribution to variance analysis is then performed on the model, but only the chosen parameters are varied, and all other parameters are held constant. There is no general rule that determines the number of parameters to include after using the MoEE, as the number of important parameters varies depending on the model. In the absence of specific information about the model, a simple cutoff could be applied, such as keeping either one percent of the total parameters, or all parameters whose total sensitivity is at least one percent of the highest global sensitivity, whichever number of parameters is greater. Arbitrary cutoffs cannot substitute for practitioner domain knowledge or experience, however, and practitioner judgment is the preferred selection method. In C

h = CF·B ·(I − A)−1f

(2)

⎧ Ri ⎪ i≤n Bi = ⎨ Pi ⎪ ⎩0 i>n

(3)

⎧ 0 j≤n ⎪ Pj ⎪ i ≤ n, j > n ⎪ 2n A i , j = ⎨ ∑m = n + 1 ∑k Ti ,m, k ·CCk + Pj ⎪ ∑k Ti , j , k ·CCk ⎪ i > n, j > n ⎪ 2n ⎩ ∑m = n + 1 ∑k Ti , m , k ·CCk + Pj

(4)

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

iteration, as the uncertainty in this factor depends only on the definition of each classification, and should not vary in space or time. Some trade data has considerable variability from year to year. We chose not to model the uncertainty of the traded mass flows of cocoa commodities, but instead to examine the relative market shares of different trade flows, normalized by cocoa content. By looking at market shares, we implicitly remove variation from year to year in country market size, and instead can base our estimates of uncertainty on the temporal variation in market share. Country market shares are calculated by dividing each trade or production flow by the total amount of cocoa-equivalent imported and grown in that year. Uncertainty was calculated from the standard deviation of the relative shares from year to year in a five-year window, and was modeled using a strictly positive normal distribution. The SI has more details and an example calculation. Characterization Factor Calculations. In LCA, we traditionally distinguish land occupation impacts, occurring during the actual use of land, and land transformation impacts, occurring due to land use change, such as clearing forests before establishment of a cocoa plantation.30 For the purpose of illustrating a new approach of sensitivity testing, we only modeled land occupation impacts in this case study, and focused on biodiversity. To assess the land occupation impact of cocoa production on biodiversity, the species richness found in cocoa plantations, Scocoa, was compared to the species richness found in nearby natural reference ecosystems, Sref. Characterization factors for land use, CF, were calculated as the relative change in species richness compared to the natural reference, as shown in eq 5:

where n is the number of countries or trade regions, i and j and row and column indices, respectively, k is trade classification indices (14 in total), h is the regionalized LCIA score vector (PDF·m2·year), f is the final demand vector (kg), CF is the vector of characterization factors of length 2n (PDF), B is the vector of computed land use flows of length 2n (m2·year/kg), R is the vector of cultivated cocoa area in countries of length n (m2), P is the vector of country cocoa production of length n (kg/year), I is the identity matrix, A is the technosphere matrix of size 2n by 2n (unitless), Ti,j,k is the mass of trade classification k traded from country i to country k (kg), CCk is the fractional cocoa content of trade classification k (unitless). Inventory Data Preparation. Data on cocoa production amounts and land areas under cultivation were drawn from the FAOSTAT database.22 Uncertainty for production amounts was classified depending on whether official statistics or FAO estimates were used and modeled using a strictly positive normal distribution, with values always greater than zero. Uncertainties for country-scale cocoa production are relatively low, as marketing boards and global agribusinesses collect highquality data.23,24 Sample standard deviations were therefore calculated using expert judgment, by multiplying the production amount by either 0.025 for official statistics or 0.05 for FAO estimates. A simple sensitivity check showed that model results did not change if these estimates were increased by a factor of 5. The normal distribution was used in the case study because it is symmetric distribution typically used in LCA and because there was not enough data to justify choosing a different distribution. The SI has more details and an example calculation. Area under cultivation was also modeled using a strictly positive normal distribution, but the standard deviation was calculated from a seven-year window over the year of interest. Because cocoa is harvested from a tree that takes five years to bear fruit,25 high variation in cultivated area from year to year is a good measure of data uncertainty. The SI has more details and an example calculation. Data on the cocoa trade came from the UN comtrade database.26 We used physical trade data in tonnes, not economic trade data in dollars, as prices vary considerably across countries and cocoa products. We used data on imports, usually collected by national customs agencies, as they are believed to have the highest available data quality.27 There are 14 trade classifications from the Harmonized Commodity Description and Coding System (HS) used for cocoa in the COMTRADE database. Some classifications, such as “Food Preparations Containing Cocoa (In Block, Slab, Bar Form; Filled)”, contain ingredients other than cocoa. We used the Codex Alimentarius28,29 to estimate the total amount and range of cocoa butter and solids in each trade classification. The mass of cocoa in each trade classification was modeled using the triangular distribution, and is detailed in the SI. In our casestudy model, the relative mass of cocoa in each trade classification is sampled once for each model Monte Carlo

CF =

Sref − Scocoa Sref

(5)

Data on the species richness of cocoa plantations and natural reference ecosystems was drawn from a literature search of peer-reviewed publications and grouped into six world regions. A total of 28 studies that compared species richness of different species groups (such as vascular plants, vertebrates, and invertebrates) for cocoa plantations and undisturbed ecosystems within one region were used. In addition, data from a global assessment of land use impacts in life cycle assessment were included.31 This study is based on data from the GLOBIO3 database,32 which consists of data from peerreviewed publications comparing the species richness and composition of different land use types with a natural reference across different world regions. We used all available cocoaspecific data, but as data availability for all six world regions was limited, we also included studies on land use types similar to cocoa such as coffee and agroforestry. For each study, we calculated a characterization factor based on eq 5. CF uncertainty distributions for each region were calculated using a weighted Gaussian kernel density estimator using Scott’s rule for bandwidth calculation.33 Data specific to cocoa was given a weight of 1, data for coffee had a weight of 0.5, and data for other agroforestry had a weight of 0.25. Weights were based on the judgment of the authors. The SI shows the calculation of the characterization factor for West Africa. The CFs represent the potentially disappeared fraction of species (PDF) in the cocoa plantations. A CF value greater than zero represents a negative impact on biodiversity, whereas D

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

Figure 1. Regionalized characterization factors (CF) for land occupation of cocoa production. The unit is potentially disappeared fraction of species (PDF). In the case study, different values are drawn for each country in each region. The highlighted points show the four values from each CF distribution used for MoEE analysis. Southeastern Asia includes Papua New Guinea.

a CF less than zero (i.e., when Scocoa > Sref) is considered a benefit to biodiversity. However, these beneficial impacts must be interpreted with caution, as the CFs just consider relative changes in species numbers, but do not report changes in the composition of species, such as a change from endemic, highly specialized forest-species to widespread generalist species. For the comparison of CFs across world regions, we need to keep in mind that the CFs represent relative changes, and therefore a 50% local loss of species in species-rich and species-poor ecosystems get the same CF, although the potential total local loss of species is much higher in a species-rich ecosystem. For a further discussion of the CFs see de Baan et al.31

of Monte Carlo model iterations is approximately 100 times the number of parameters assessed. Case Study Results. Figure 1 shows the regionalized characterization factors for land use due to cocoa production. Considerable uncertainty is present in all regions. The Kruskal−Wallis test indicated that the CF distributions are independent (p ≪ 0.001), and did pairwise Kolmogorov− Smirnov tests for each possible pairing (p < 0.001). Negative CF values occur when more species are recorded on cocoa plots than in the reference state. Although data from coffee and other agroforestry was used in addition to data specific to cocoa cultivation, none of the extremes in any region are due to the inclusion of the non cocoa-specific data. The SI data file has all the data used to generate these characterization factor distributions. Figure 2 shows results from the MoEE analysis. Characterization factors had the highest global sensitivity, and there is a clear difference between the ratio of absolute average to standard deviation for CF parameters and almost all other parameters. The points in Figure 1 show the values used in MoEE analysis, and in five of the six regions considered, at least one of the representative points is negative. Increases in trade from or production in a certain country could therefore lead to decreases in the model score if the CF in that country sampled during Monte Carlo iteration was negative. The changes in model scores caused by changes in cultivated area, yield, or trade flows could therefore be positive or negative, depending on the sign of the CF, and the standard deviation of the model score changes was higher for these parameters than for the CF parameters.



RESULTS Convergence of CTV Scores. We used results from the CTV numerical simulation to fit eq 6, which relates the average squared ROCC, uncorr, for a parameter uncorrelated with the model outputs to the number of Monte Carlo model iterations, n: uncorr =

1.042 n1.004

(6)

Equation 6 cannot be used directly in CTV calculations, but rather provides an estimate of the expected CTV error for parameters not correlated with model results. Because uncorr is never zero, if there are a large number of uncorrelated parameters, each with a squared ROCC of uncorr, the CTV of an important parameter will be small if n is relatively small. SI Figure S2 shows CTV scores converge only when the number E

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

Indonesian CF’s CTV score is due to the 80% increase in market share of Indonesian cocoa in Japan from 2000 to 2005. Although there were only few large changes in the CTV results of Switzerland or Japan from 2000 to 2005, there was a marked increase in the importance of Ivory Coast production to model output variability in the U.S. There were small increases in the proportion of Ivory Coast cocoa, and small decreases in proportion of Indonesian cocoa, in the U.S. market. However, the change in CTV results was due to almost 50% larger and more uncertain direct trade flows from the Ivory Coast to the U.S. This increase in direct importation is consistent with the entry of U.S.-based agribusinesses such as ADM into West Africa24 during this time. Direct trade flows cause more model variance than independent indirect trade flows, whose variability is significantly reduced through the averaging effect.



DISCUSSION Two-Step Sensitivity Analysis. The two-step sensitivity analysis proposed in this paper has a number of advantages over CTV alone. First, it can lead to a significant reduction in the total number of model evaluations. The precise reduction is a factor of how many MoEE trajectories are used, and how precise the desired CTV values should be, but a reduction of an order of magnitude or more is expected. This reduction comes about because the MoEE normally needs 10 model evaluations per parameter, as 10 trajectories are computed, and CTV only converges with approximately 100 model evaluations per parameter. The number of model evaluations needed for subsequent CTV analysis with the reduced parameter set is insignificant. Using both the MoEE and CTV also gives two separate variance-based measures of parameter sensitivity. The strengths of both methods are complementary. The MoEE uses trajectories that attempt to fill the parameter hypercube, because they are chosen by maximizing their distances from each other. This means that the MoEE, like Latin Hypercubic sampling, is good at evaluating models for all possible combinations of the selected parameter values, and especially for combinations of extreme parameter values. However, the MoEE is explicitly a screening method, and should be followed with a more detailed sensitivity analysis for the selected parameters. CTV draws random samples from each probability distribution, not from a restricted set of possible values, and ensures more complete coverage of all parameter values and value combinations. CTV results are also easier to interpret and communicate. The MoEE assesses both parameter sensitivity and nonlinearity, which is very helpful when trying to interpret and improve complex parametrized LCA systems. Manual sensitivity analysis or algorithms which assume strictly linear systems cannot handle the hundreds or thousands of interacting parameters which will be in the next generation of inventory databases such as ecoinvent 3. Because the MoEE examines all parameters and complete parameter ranges, it can also help detect errors in parametrized inventories that would not be otherwise observed. The MoEE can also be applied to groups of similar parameters. Examining groups of parameters can reduce the number of model evaluations needed, which would be helpful for computationally expensive LCA models. However, the use of groups instead of individual parameters can make it more

Figure 2. Method of elementary effects results plot for market mix of 1 kg of cocoa-equivalent in the U.S. in the year 2000, and shown in larger, bold circles. Results for Japan and Switzerland in 2000 and 2005 and the U.S. in 2005 are also shown, but with smaller markers.

The 100 model parameters with the highest μ* scores were selected for further analysis with CTV. Table 1 shows the CTV Table 1. Contribution to Variance (CTV) Scores for Top Three Parameters for One Kilogram of Cocoa-Equivalent Consumed from the Market in Switzerland, the U.S., and Japan in 2000 and 2005 Switzerland parameter

U.S.

Japan

CTV

parameter

CTV

parameter

CTV

Ghana CF Ivory Coast CF Ecuador CF

0.44 0.29

0.42 0.31

Ghana CF Indonesia CF

0.52 0.18

0.1

2000 Indonesia CF Ivory Coast CF Ghana CF

0.11

Ivory Coast CF

0.12

Ivory Coast CF Ghana CF Ecuador CF

0.39

2005 Ivory Cost CF

0.64

Ghana CF

0.44

0.38 0.095

Indonesia CF Ghana CF

0.16 0.053

Indonesia CF Ivory Coast CF

0.32 0.12

scores for the top five parameters for the consumption of one kilogram of cocoa-equivalent from the market mixes of Switzerland, the U.S., and Japan in 2000 and 2005. These were the first and last years that could be used with the uncertainty modeling, which requires three years of data before and after the chosen year for time series analysis. In five of the six cases, the characterization factors for land use impacts in the Ivory Coast and Ghana contributed to more than 50% of the observed model variance. In Switzerland, the relative decrease in the CTV of the Ghanaian CF from 2000 to 2005 is due to a decline in the relative share of Ghanaian production in the Swiss market. In Japan, the market shares of Ivory Coast and Ghana are almost equal and do not change significantly throughout time, but the combination of more uncertain trade and agricultural data from Ghana lead to its characterization factor having a higher CTV score. The increase in the F

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

Implications. Sensitivity analysis will face a trade-off between comprehensiveness and computation speed. Although regionalization will increase the size of inventory databases, regionalized LCA calculations with complicated geospatial data can still be done using a modification of standard LCA matrix math.9 However, parametrized data sets will be much slower to compute, as the speed advantages of tuned linear algebra libraries will be lost. An order of magnitude increase in calculation speed may not be enough if the regionalized and parametrized inventories are multiple orders of magnitude slower to evaluate. One possibility for future sensitivity analysis is a tiered system: inventory database developers use the MoEE to understand the behavior of and interaction between parametrized data sets, and to calculate in advance the list of sensitive parameters for different impact assessment methods. These developers could also use CTV to identify sensitive parameters, and to allocate resources to reduce overall variance in LCA results. LCA practitioners could then perform twotiered sensitivity analysis against a reduced set of foreground and selected background parameters in a reasonable amount of computation time. Software developers can use the open source software developed for this article either directly or as a starting point for the integration of the MoEE and CTV into software tools. Life cycle inventories are starting to provide support for parametrization, including dependent equations with multiple parameters, and regionalized data sets. Regionalized impact assessment methods are also increasingly being developed. These advances allow for more powerful LCA models, but also increase the need for sensitivity analysis techniques that can help assess these complicated systems. The two-step approach described in this manuscript allows for sensitivity analysis of LCA models with many parameters, and does not make assumptions about parameter distributions or interrelations. By providing information on parameter global sensitivity, nonlinearity, and contribution to overall model variance, this new approach can help interpret and improve complex LCA system models.

difficult to identify all the important parameters because of averaging effects.19 De Koning et al20 applied the Fourier amplitude sensitivity test (FAST), another variance-based sensitivity test. FAST could be used instead of CTV, and is a popular and wellunderstood way to calculate sensitivity indices. CTV is much simpler to calculate and implement, however, and provides results that are easier to interpret. Further comparative discussion of global sensitivity analysis techniques can be found in Saltelli et al.1 Case Study. In the case study, the land use characterization factors were the dominant uncertain parameters. In all but one of the countries and years examined, the characterization factors for land use due to cocoa production in Ivory Coast and Ghana accounted for over 50% of the total model variability. This variability came from the wide range of the empirical characterization factor distributions, most of which even crossed the origin, and therefore could be positive or negative. This variation is a result of several different factors. In the case study, we combined data from different types of cocoa production systems, such as high-intensive full-sun and lowintensive agroforestry systems. These different production systems have been shown to have diverging impacts on biodiversity.34 The calculated CFs also combined data of different species groups, and aggregated local effects to large regions covering different biomes. To reduce uncertainty of the case study, production system-specific inventory data and characterization factors for different cocoa production systems and biomes would be required. Despite the uncertainties in trade and country agricultural data, the use of market mixes and cocoa-content normalization made models of the trade of cocoa feasible. Although the trade model is simplified, it can be used in cases where information on domestic stocks and consumption and country-specific value chains for the transformation of trade classifications are not available. The case study also shows the importance of considering uncertainty in impact assessment, which is rarely done and not supported in most LCA software. Uncertainty in the CFs was the dominant driver of model output variability in the case study. Regionalized LCA will necessitate software support for different CF probability distributions being applied to many different regions. Uncorrelated error in CTV. The uncorrelated error term in CTV derives from the use of the square of the rank-order correlation coefficient. Although the ROCC of two uncorrelated data sets is normally distributed and has an expected value of zero,16 the squared values will always be positive, and therefore the expected value of the squared ROCC of two uncorrelated data sets is always positive. The CTV equation requires that the ROCCs be squared so that each numerator is positive and the denominator does not include the summation of negative numbers. Other nonparametric correlation coefficients, such as a weighted rank-order correlation coefficient35 or Kendall’s Tau,36 could be used to calculate CTV scores. Both of these correlation coefficients are reported to have higher discriminating power than Spearman’s correlation coefficient. They would not eliminate the uncorrelated error term, however, because the CTV equation would still require squaring their values, leading to a nonzero expected correlation for uncorrelated data.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information has details of the computer code developed for this manuscript, discussion of the mathematical properties of contribution to variance, the transformation of uncertainty distributions to the unit interval, detailed case study calculations and CF data, and a graphical illustration of the mathematical model. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*Phone: +41-44-633-71-45; fax: +41-44-633-10-61; e-mail: [email protected]. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest. G

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology



Article

(21) Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S. Global Sensitivity Analysis: The Primer; John Wiley & Sons: Chichester, England, 2008. (22) FAOSTAT Statistical Database. Food & Agriculture Organization: Rome, Italy, 2012; http://faostat.fao.org/. (23) Dand, R. the International Cocoa Trade, 3rd ed; Woodhead: Oxford, England, 2011. (24) Musselli, I., Combe, O. Cocoa Study: Industry Structures and Competition; UNCTAD: New York, NY, 2008; http://r0.unctad.org/ infocomm/Studies/Cocoa_Study_2008.htm. (25) Willison, K. Coffee, Cocoa and Tea; CABI Publishing: Wallingford, England, 1999. (26) UN comtrade. United Nations: New York, NY, 2012; http:// comtrade.un.org/. (27) Bouwmeester, M., Oosterhaven, J. Technical Report: Inventory of Trade Data and Options for Creating Linkages. Report of the EXIOPOL Project; University of Groningen: Groningen, Netherlands, 2007; http://www.feem-project.net/exiopol/userfiles/file/Cluster%20III_ M18_PU/EXIOPOL_DIII_1_a_3_final.pdf. (28) Draft Revised Standard for Cocoa Powders (Cocoas) and Dry Mixtures Of Cocoa and Sugars; Food & Agriculture Organization: Rome, Italy, 2001. (29) Codex Standard for Chocolate and Chocolate Products; Food & Agriculture Organization: Rome, Italy, 2003. (30) Milà i Canals, L.; Bauer, C.; Depestele, J.; Dubreuil, A.; Knuchel, R.; Gaillaird, G.; Michelsen, O.; Müller-Wenk, R.; Rydgren, B. Key elements in a framework for land use impact assessment within LCA. Int. J. Life Cycle Assess. 2007, 12 (1), 5−15. (31) de Baan, L.; Alkemade, R.; Koellner, T. Land use impacts on biodiversity in LCA: A global approach. Int. J. Life Cycle Assess. 2012, DOI: 10.1007/s11367-012-0412-0. (32) Alkemade, R.; van Oorschot, M.; Miles, L.; Nellemann, C.; Bakkenes, M.; Ten Brink, B. GLOBIO3: A framework to investigate options for reducing global terrestrial biodiversity loss. Ecosystems 2009, 12 (3), 374−390. (33) Scott, D. Multivariate Density Estimation: Theory, Practice, and Visualization; John Wiley & Sons: Chicester, England, 1992. (34) Donald, P. Biodiversity impacts of some agricultural commodity production systems. Conserv. Biol. 2004, 18 (1), 17−37. (35) da Costa, J.; Roque, L. Limit distribution for the weighted rank correlation coefficient. REVSTAT − Stat. J. 2006, 4 (3), 189−200. (36) Kendall, M. A new measure of rank correlation. Biometrika 1938, 30 (1/2), 81−93.

ACKNOWLEDGMENTS We deeply thank the three anonymous reviewers whose detailed comments greatly improved this manuscript.



REFERENCES

(1) Saltelli, A., Tarantola, S., Campologno, F., Ratto, M. Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models; John Wiley & Sons: Chicester, England, 2004. (2) Bjorklund, A. Survey of approaches to improve reliability in LCA. Int. J. Life Cycle Assess. 2002, 7 (2), 64−72. (3) Huijbregts, M.; Norris, G.; Bretz, R.; Ciroth, A.; Maurice, B.; von Bahr, B.; Weidema, B.; de Beaufort, A. Framework for modelling data uncertainty in life cycle inventories. Int. J. Life Cycle Assess. 2001, 6 (3), 127−132. (4) ISO. 14044:2006 - Environmental ManagementLife Cycle AssessmentRequirements and Guidelines; ISO: Geneva, Switzerland, 2006. (5) Saltelli, A.; Annoni, P. How to avoid a perfunctory sensitivity analysis. Environ. Model Software 2010, 25 (12), 1508−1517. (6) Heijungs, R. Identification of key issues for further investigation in improving the reliability of life cycle assessments. J. Clean. Prod. 1996, 4 (3−4), 159−166. (7) Heijungs, R.; Kleijn, R. Numerical approaches towards life cycle interpretation: Five examples. Int. J. Life Cycle Assess. 2001, 6 (3), 141− 148. (8) Sakai, S.; Yokoyama, K. Formulation of sensitivity analysis in life cycle assessment using a perturbation method. Clean Technol. Environ. Policy 2002, 4 (2), 72−78. (9) Mutel, C.; Pfister, S.; Hellweg, S. GIS-based regionalized life cycle assessment: How big is small enough? methodology and case study of electricity generation. Environ. Sci. Technol. 2012, 46 (2), 1096−1103. (10) Heijungs, R. Sensitivity coefficients for matrix-based LCA. Int. J. Life Cycle Assess. 2010, 15 (5), 511−520. (11) Geisler, G.; Hellweg, S.; Hungerbuhler, K. Uncertainty analysis in life cycle assessment (LCA): Case study on plant-protection products and implications for decision making. Int. J. Life Cycle Assess. 2005, 10, 184−192. (12) Milà i Canals, L.; Azapagic, A.; Doka, G.; Jefferies, D.; King, H.; Mutel, C.; Nemecek, T.; Roches, A.; Sim, S.; Stichnothe, H.; Thoma, G.; Williams, A. Approaches for addressing life cycle assessment data gaps for bio-based products. J. Ind, Ecol. 2011, 15 (5), 707−725. (13) Cooper, J.; Noon, M.; Kahn, E. Parameterization in life cycle assessment inventory data: Review of current use and the representation of uncertainty. Int. J. Life Cycle Assess. 2012, 17 (6), 689−695. (14) Weidema, B., Bauer, C., Hischier, R., Mutel, C., Nemecek, T., Vadenbo, C., Wernet, G. Overview and Methodology: Data Quality Guidelines for the Ecoinvent Database Version 3. Final Draft Revision 2; Ecoinvent Centre: Dübendorf, Switzerland, 2011; http://www. ecoinvent.org/fileadmin/documents/en/ecoinvent_v3_elements/01_ DataQualityGuideline_FinalDraft_rev2.pdf. (15) Pfister, S.; Koehler, A.; Hellweg, S. Assessing the environmental impacts of freshwater consumption in LCA. Environ. Sci. Technol. 2009, 43 (11), 4098−4104. (16) Lyerly, S. The average spearman rank correlation coefficient. Psychometrika 1952, 17 (4), 421−428. (17) Cureton, E. The average spearman rank criterion correlation when ties are present. Psychometrika 1958, 23 (3), 271−272. (18) Morris, M. Factorial sampling plans for preliminary computational experiments. Technometrics 1991, 33 (2), 161−174. (19) Campolongo, F.; Cariboni, J.; Saltelli, A. An effective screening design for sensitivity analysis of large models. Environ. Model. Software 2007, 22 (10), 1509−1518. (20) de Koning, A.; Schowanek, D.; Dewaele, J.; Weisbrod, A.; Guineé, J. Uncertainties in a carbon footprint model for detergents; quantifying the confidence in a comparative result. Int. J. Life Cycle Assess. 2010, 15, 79−89. H

dx.doi.org/10.1021/es3050949 | Environ. Sci. Technol. XXXX, XXX, XXX−XXX