Scenario analysis and sensitivity exploration of the MEDEAS Europe energy-economy-environment model

decision-makers rely heavily on Integrated Assessment Models to guide the decarbonisation of the en- ergy system. Uncertainty is embedded in the assumptions these models are built upon. Unless those uncertainties are adequately assessed, using Integrated Assessment Models for policy design is unadvised. In this work we run Monte Carlo simulations with the MEDEAS model at European Union scale to assess how the uncertainties on the main drivers of the transition affect key socioeconomic and environmental indicators. In addition, One-at-a-time sensitivity exploration is performed to grade the contribution of a set of model parameters to the uncertainty in the same key indicators. The combination of the uncertainties in the model drivers magnify the uncertainty in the model outputs, which widens over time. Parameters affecting sectorial and households ’ energy ef(cid:51)ciency and households ’ transport energy use ranked amongst the most impacting ones on simulation results.


Introduction
In the process to achieving long-term sustainability of human societies, and with capitalism as the global economic system of choice, all hopes are put on the urgent replacement of fossil fuels by renewable energy as the mean to avoid the worst impacts of the current climate crisis and the foreseeable energy crisis.
Compared to past global challenges such as the ozone layer depletion, maintaining the status quo while transitioning our societies in such context of climate change, environmental degradation and potential energy scarcity is undoubtedly amongst the most complex humanity as a whole has ever faced.
While certain critical planetary boundaries have already been crossed [1,2], putting the climatic stability and resource abundance that our societies have enjoyed during the Holocene at stake, failing to take the right decisions in the coming years may bring our planet and our civilization to a point of no return. In such challenging and uncertain context, all kinds of tools, perspectives and approaches will be needed to guide decision-making. In this sense, numerical models, and particularly Integrated Assessment Models (IAMs) are a group of mathematical models used to portray the social, economic, environmental, climatic and institutional dimensions of, in the context of this work, the energy transition [3]. As such, they represent valuable decision support tools, since they provide a riskless manner of exploring alternative scenarios and policies [4,5].
Despite the great advances made in the development of IAMs in the last decades, most of them still share a core set of common assumptions whose validity is being disputed in scienti c forums. First, IAMs are generally characterized by a rather sequential structure with limited feedbacks among the represented subsystems. It is especially relevant the omission of climate change impacts [6][7][8]. Second, a lack of plurality in the methods to represent the economic dimension has been detected in the literature, dominated by assumptions of conventional equilibrium through optimization methods, aggregated production functions as well as the widespread use of prices as indicators of scarcity [9]. Third, fossil fuel resource abundance, understood as the vast geological availability accessible at an affordable price, is a default assumption in most of the prominent IAMs used for climate policy analysis. Also, it is usually assumed that the resource base of renewable energy sources (RES) provides no practical limitation if adequate investments are forthcoming [10]. Finally, most models disregard the implications that the future energy and mineral investments to achieve the transition to renewables may have for the system [11].
The MEDEAS models were developed to ll those gaps. Leaving technological and economic parameter differences and licensing aside, based on the classi cation used in Ref. [8], the MEDEAS models are comparable to the ETM [12], LEAP [13], NEMS [14], POLES [15], RETScreen [16] and IEA's WEM [17] models. However, the following aspects set the MEDEAS models apart from their counterparts. First, their implementation using the system dynamics approach allows to represent none-equilibrium processes and feedback loops (e.g. climate damage); second, biophysical limits in terms of energy and mineral resources are enforced: the availability of fossil fuels is represented using Hubbert curves [18], and the minerals consumption is referred to their actual resources and reserves; third, by using physical variables (energy) in the economic module, potential energy scarcity events are represented without losing sight of the economy (by using Input Output Tables and energy intensities [19]); and nally, the Python version of the MEDEAS models (pymedeas) is distributed for free and is fully open source.
Despite the unquestionable potential of IAMs, the largely uncertain factors that will shape the future energy system, including technological innovation, resource availability, socio-economic dynamics and geopolitics [20], are unavoidably embedded in the thousands of hypothesis and parameters included in these models and propagate to their outputs [21]. Consequently, unless suf cient testing of their projections is provided, these models risk losing credibility [22,23]. Most importantly, scienti c models are today's intermediaries between science and policy [24], hence uncertainty represents a non-negligible concern [25] that has to be addressed in order for the models to serve society [26].
Uncertainty and Sensitivity analysis are methodologies that tackle that particular issue. Uncertainty Analysis (UA) aims at assessing the uncertainty of the model projections resulting from the uncertainty in the model inputs, not trying to identify its origin, while Sensitivity Analysis (SA) is used to identify the input parameters' relative contribution to the outputs uncertainty. To this regard, SA is not to be considered as an alternative to UA but rather as its complement [27].
The need to include UA and SA in the modelling process is backed by a large number of publications that count these techniques among the best modelling practices and especially for their application in policy design [20,[28][29][30]. The information extracted from their application provides a better understanding of the functioning and behavioural boundaries of the model. Most importantly, these techniques must make it possible to explain results from the model in terms of the mechanisms that drives the model [30].
Although the application of UA and SA in energy-economy IAMs is limited [24,30], examples include: Ref. [22], who performed Global Sensitivity Analysis (GSA) of Res-IRF, an energy-economy model of the demand for space heating in French dwellings; Ref. [31] used UK energy system model, ESME, to explore trade-offs in cost effective energy transition scenarios and performed a sensitivity analysis to explore the uncertainties that have most impact on the transition; Ref. [32] performed a sensitivity analysis on the GET 7.0 model to analyse how the development of the energy supply system in a carbon-constrained world in uences the cost-effectiveness of fuels and propulsion technologies in the transportation sector. Other examples of the application of SA techniques in IAMs are reviewed in Ref. [33]. Based on the low number of such works, and taking into account the vast and largely unknown number of parameters that IAMs include, the knowledge gaps and uncertainties on the shape the energy transition will take remain even larger.
To date no study has been published dealing with uncertainty in the MEDEAS models. Only in Ref. [34] a scenario analysis was performed, though in the form of a parametric optimization. Hence, the current work represents the rst application of uncertainty and sensitivity procedures to the MEDEAS models, and a further claim for the need to apply these techniques in IAMs used in policy design.
Additionally, in this work parameter sensitivity is analysed both qualitatively and quantitatively, and 3 different measures of sensitivity are combined to derive results. The rst measure consists on the use of Spider plots, which provides visual clues on the impact of each input on each output. On the other hand, the quantitative measures are the RMSD and the Euclidian distance. The rst provides an indication of the sensitivity of each output parameter relative to each input along the analysed timeframe, while the second measures the impact of each input on the combined trajectory of all the selected outputs altogether.
By using these procedures, the main objectives of this work are (1) to assess the uncertainties embedded in the model projections and (2) to identify the main drivers of the energy transition of the EU28, hence provide clues on the key aspects that need to be dealt with in order to make the energy transition happen. After the Introduction, this document is structured as follows: in the Methods section we rst present the theoretical background of the MEDEAS models. After that we describe the simulations that were run and present how the results were processed in order to obtain the indicators used to assess uncertainty and sensitivity. The Methods section is followed by the Results, the Discussion and the nal conclusions.

The MEDEAS models
MEDEAS (Modelling sustainable Energy system Development under Environmental And Socioeconomic constraints) is a Horizon 2020 research project that started in 2015 with the aim of developing new IAMs to tackle the aforementioned limitations of most mainstream IAMs. This project represents a joint effort of a consortium of 12 European institutions, including public and private research centres and universities, energy agencies and IT services companies, to provide policy makers with a modelling tool to test new and existing policies to achieve a more sustainable European Energy system.
More speci cally, the MEDEAS are a set of a policy-simulation dynamic-recursive models sharing the same conceptual modelling approach and designed applying System Dynamics. They are available at three geographical scales: World, EU28 and country-level (Austria and Bulgaria). The one used in this work is that at European scale.
Although the models were originally developed in the proprietary Vensim® software, the of cial version of the MEDEAS models are written in Python and distributed free of charge under the MIT license, guaranteeing transparency, reproducibility and traceability [4].
All three models are structured in seven main conceptual submodules: Economy, Energy availability, Energy infrastructures, Materials, Land-use, Climate/Emissions, and Social & Environmental impact indicators ( Fig. 1).
The MEDEAS models dynamically operate as follows: for each period, a sectoral economic demand is estimated from exogenous pathways of expected Gross Domestic Product per capita (GDPpc) and population evolution. The nal energy demand required to ful l production is obtained using energy-economy hybrid input-output analysis, and energy intensities by type of nal energy. The energy sub-module computes the available nal energy supply, which may or may not satisfy demand, adapting the economic production to the available energy. The materials required by the economy, with emphasis on those required by alternative green technologies, are estimated; this allows to assess eventual future mineral bottlenecks. The new energy infrastructure requires energy investments, whose computation allows to estimate the variation of the EROI (Energy Return over Energy Invested) of the system, which in turn affects the nal energy demand. The climate submodule computes the greenhouse gas (GHG) emissions associated to the resulting energy mix (complemented by exogenous pathways for nonenergy emissions), which feeds back to the economy, affecting nal demand. Additional land requirements are accounted for. Finally, the social and environmental impacts are computed. For more detail the reader is referred to Refs. [36,37].
Despite their short existence, the MEDEAS models have already been used to study different aspects of the energy-economy-environment interrelations in 6 publications. In Ref. [36] the MEDEAS model at World scale was presented, and the main features and hypothesis were discussed. In Ref. [37], the Python version of the EU28 model (pymedea-s_eu) was introduced, the principles behind its development approach were laid-out (openness, transparency, user friendliness and community-based design) and the effect of the main features of the model (evolution of sectorial energy ef ciencies, EROI (Energy Return on Investment) feedback, climate change impacts and fossil fuels availability) was demonstrated through a scenario analysis. In the third study [38], the MEDEAS model at World scale was used to evaluate the energy (based on the EROI) and the material costs of the transition. Ref. [34] performed an optimization study aimed at nding the values of a set of model parameters in order to t the simulated CO 2 emissions to previously obtained decarbonisation pathways that would allow to remain below 2 • C of global warming [39]. The fth work [19] presents the macro-economic module included in the MEDEAS models and uses simulation results obtained from different input scenarios to highlight the con ict between economic growth, climate policy and the sustainability of resources. Finally, Ref. [40] describes a novel methodology, implemented in the MEDEAS models, to estimate energy demand based on the projection of the evolution of sectoral nal energy intensities.

Simulations
The simulations are divided between those used for the Uncertainty Fig. 1. Graphical representation of the modules of the MEDEAS models [35]. Analysis (UA) and those used for the Sensitivity analysis (SA). Thus, in the rst group of simulations (UA), the propagation of uncertainty from inputs to outputs is assessed, while in the second group (SA), the contributions of the inputs to the uncertainty of each of the outputs is analysed.
The protocol presented herein corresponds to running multiple simulations with different combinations of the input parameters values, each within a prede ned range. This way, uncertainty in the inputs is propagated to the outputs and can be assesses/quanti ed.
All simulations comprise the time-frame between 1995 and 2050, with a time-step of 0.03125 years (11.4 days). The outputs are analysed at an annual frequency. Model projections start from 2015, while before that the presented results correspond to historical data.
The scenario used for all simulations is the Business as Usual (BAU), which perpetuates historical trends of economic growth, sustained by fossil fuels, in which little efforts are put into making the transition to decarbonizing the energy system. For further details, the description of the most relevant inputs and assumptions that characterize the BAU narrative in the MEDEAS models can be found in Ref. [40] (Appendix C and Table C. 1).

Uncertainty Analysis
The 27 input parameters perturbed during the scenario analysis belong to a subgroup of the exogenous model parameters that the user is expected to tweak to create different simulation scenarios (model drivers). Their names, description and their assigned minimum and maximum values are presented in Table 1.
Among all parameters that de ne scenarios in the MEDEAS models, these 27 were selected due to their uncertain nature and because they were expected to have a large impact on the simulation outputs. The list includes the desired GDPpc annual growth rate, the annual population change, and the techno-ecological potentials and growth of the deployment rates of the different RE technologies included in the model (solar photovoltaic (PV), Concentrated Solar Power (CSP), onshore and offshore wind, etc.) (see Table 1). The upper and lower values of the ranges for each of them were selected making a trade-off between covering the widest possible range of the values to check the model stability, and limiting the values to realistic gures for the uncertainty analysis.
Version 1.2 of the Vensim® implementation of the model (MEDEA-S_EU v1.2) was used to take advantage of the built-in Sensitivity Simulations functionality, which in this work was used to evaluate uncertainty. One thousand multivariate Monte Carlo simulations were run on a conventional workstation, sampling from random uniform distributions between the minimum and maximum values given to each input parameter (Table 1).

Sensitivity exploration
For the sensitivity exploration, 19 xed-value (mostly based on literature values) model parameters were selected. The focus was put, again, on those parameters whose values were judged more uncertain. The minimum and maximum values they were given were based on Table 1 Input parameter names (Vensim® nomenclature), description and minimum and maximum values used to obtain the random uniform distribution for each of them for the scenario analysis. expert best guesses and are also presented in Table 2.
In addition, to evaluate the impact of two endogenous parameters on the model outputs, they were both multiplied by a constant, that either divided (minimum value of the range) or multiplied (maximum value) the actual value of the variable by 2. The two endogenous variables were in fact vectors, and only some of their dimensions were perturbed (see Table 3).
The Sensitivity exploration was carried out using version 0.3.0 of pymedeas_eu model, which is the Python translation of the Vensim® model used for scenario analysis. Forty-three simulations were run, modifying one parameter at a time (OAT): 21 with the upper values of the range of each parameter, 21 with the lower values and 1 with nominal values (see Tables 2 and 3).

Sensitivity indexes and results processing
The list of endogenous parameters (model outputs) upon which the uncertainty of the inputs is assessed, was confectioned based on a qualitative comparison of the MEDEAS model with two of the most widely used energy-economy models (TIMES [41] and LEAP [13]), and includes the parameters that are present in all three models (Table 4). This previous work was done with the aim of facilitating a future comparison between the outputs of the three models at country scale (Austria and Bulgaria).

Uncertainty analysis
Results of the scenario analysis were qualitatively interpreted using fan charts showing the con dence bounds around the mean of the distribution of each output variable for the 1000 Monte Carlo simulations.

Sensitivity exploration
On the other hand, the sensitivity measures (qualitative and quantitative) used to analyse the outputs of the 43 simulations consisted on: a) plotting the values of all output variables at the end of the simulation (year 2050) against the percentage of change of each perturbed parameter with respect to its nominal value (spider plots). b) the Root Mean Square Deviation (RMSD) between the curves obtained for each output with the minimum and maximum values of the perturbed input parameters. c) the Euclidean distance between all aggregated outputs obtained with perturbed and nominal values of each input parameter (maximumnominal and minimum-nominal distances).
The expressions below de ne the input (Eq. (1)) and output (Eq. (2)) parameter spaces, that will be used in subsequent equations: Where p i corresponds to any of the input parameters of Table 1 and z j refers to any of the output variables of Table 4 and u j are any other parameters or variables not included in the input parameter space P that each output variable depends on. Note that not all output variables depend on all input parameters in P. t is time (Eq. (3)).
The RMSD between the output variables obtained with the minimum and maximum values of each perturbed parameter (RMSD p i j ) are calculated with Eq. (4): The RMSD is divided by the percentage change of the perturbed parameter (Δp i ) so that all RMSD obtained for the same output variable (perturbing different inputs) are comparable among them (Eq. (5)). To ensure that all outputs have the same importance on the calculation of the Euclidean distance, they are previously standardised using Eq. (6): where f j is the time-series to standardise (endogenous variable), μ j is its mean value and σ j is its standard deviation. The standardised time-series (f j ) will have a mean of 0 and a standard deviation of 1. The Euclidean distance (scalar) between the vectors of normalised outputs obtained with the maximum and nominal values of the input parameter p i , at time t (E pi,max−nom t ) is obtained with Eq. (7): Where, Ž

Table 4
List of common outputs of the MEDEAS_EU and pymedeas_eu models (using the naming convention of pymedeas_eu) to be analysed in scenario analysis and sensitivity exploration. To get a single value for the whole simulated period, the Euclidean distances at all simulation times are added as follows (Eq. (8)): The resulting value is divided by the percentage of change of the perturbed parameter so that the Euclidean distances obtained for every input parameter are comparable (Eq. (9)): The same procedure is done to obtain the time aggregated Euclidean distance between the minimum and nominal values of each input parameter p i (E p i ,min−nom ), which are also divided by the percentage change of the minimum and nominal values of the input parameter (δp min−nom i ).

Uncertainty analysis
None of the 1000 Monte Carlo simulations failed to converge. Therefore, the minimum and maximum values of the input parameters shown in Table 1 correspond to the tested stability ranges of the model.  Table 4 are presented.
The values of the standard EROI obtained from the Monte Carlo Experiments start at ca. 11 in 2015 and by 2050 the spread goes from as low as ca. 7 to as high as 13.5 (Fig. 2). The distribution is uniform with a slight negative skew. On the other hand, the distributions for GDP, GDPpc, total nal energy consumption (TFEC) and the CO 2 emissions are very similar (Fig. 3) and show a very large positive skew. Some of the simulations for these 4 output variables display positive exponential trends. The share of electricity produced from RES shows a large uncertainty with Platykurtic distribution (negative kurtosis) and a negative skew towards the end of the simulated period (Fig. 4). Finally, the uncertainty of the land requirements is again very large (Fig. 5), and in this case the top limit is the result of the decreasing land availability which acts as a ceiling for further renewable capacity installation.

Sensitivity exploration
Spider plots are used to qualitatively assess the impact of the perturbation of the parameters from Tables 2 and 3 on the values of the  output variables from Table 4   Input parameter min_energy_intensity_vs_initial represents the minimum value that the energy intensity of each economic sector can potentially reach, with respect to their intensities in 2009. Perturbations above and below its nominal value (30%) also produced noticeable changes in the output variables, being the system's EROI and the total nal energy consumption (TFEC) the output variables more signi cantly impacted. CO 2 emissions also increased for higher values of this parameter.
Parameter threshold_remaining_potential_new_capacity, corresponds to the threshold value of the remaining potential capacities of the RE for electricity generation below which the planning of new renewable infrastructure starts to decline due to decreasing returns. Increasing the nominal value of this parameter (0.5) had a negative impact on the share of RES and reduced land requirements.
Input parameter a1_coef_th, with units of EJ/10 12 1995 US$, is a coef cient involved in the calculation of the variation of the energy intensities of the household transport sector resulting from using different energy carriers (gas, liquids and electric batteries), and had a nominal value of 1.46. This parameter had a negative relation with TFEC, CO 2 emissions and the standard EROI of the system. CO 2 emissions were greatly affected by the minimum value taken by the capacity factor (min_cp_nuclear, nominal value = 0.6), which represents the average fraction of the time that nuclear plants are producing energy with respect to the total time. The value of this parameter may tip the balance between requiring new nuclear capacity or not. Tables 5-7 show the highest three values of the normalised RMSD between each output variable obtained with the minimum and maximum values of the range of each input parameter. For each output variable (columns in the tables), the higher the value of the normalised RMSD, the more sensitive it is to changes in the perturbed parameter (rows in the tables). Table 8 shows to which perturbed parameter each output variable was most sensitive, based on the values of the normalised RMSD. Table 9 shows that input parameter min_energy_intensity_vs_initial was in the top three of the most impacting parameters for all 15 outputs. Variable variation_nonxdashxenergy_use and parameter a1_coef_th came second and third, being 8 and 7 times among the top three, respectively. These results con rm the qualitative results obtained with the spider plots.
Similar conclusions were drawn based on the values of the timeaggregated Euclidian distances of all outputs obtained with the minimum and nominal and maximum and nominal values of each input parameter (Fig. 7). Accordingly, the most impacting of the selected input parameters was the min_energy_intensity_vs_initial, followed by a1_conf_th and min_energy_intensity_vs_initial_h. Variable variation_nonxdashxe-nergy_use came in at the 4 th place.
On the other hand, perturbations in parameters min_lifeti-me_ev_batteries, exponent_availability_conv_gas, share_max_of_change_ vs_historical_mean_rate and share_max_of_change_vs_historical_ mean_rate_h (see descriptions in Table 2) did not affect the analysed output variables in any particular way. Other parameters that had little impact on the outputs were share_res_elec_generation_ curtailedxandxsored, share_gasxdivxcoalplusxgasx_for_heat_plants and esoi_phs_depleted_potential.

Discussion
As Ref. [42] points out, it must be noted that sensitivity and uncertainty analysis are themselves uncertain, because there is considerable uncertainty in quantifying the uncertainty in input factors. Hence, the criteria for the selection of the input parameter ranges must be transparent. In this work, the minimum and maximum values of the input  parameters for the uncertainty analysis (Table 1) were selected using two criteria: they should be suf ciently far apart to allow testing the model for stability, but at the same time remain realistic. Most of the parameters of the list correspond to the techno-ecological potentials for the different RES technologies, whose values vary orders of magnitude depending on the bibliography. The list also includes the annual increase in the installed capacity of the same technologies which, at least on the rst stages of the transition, will not be affected by biophysical (materials or energy scarcity) or technical constraints (grid integration) and hence have been assigned wide ranges. Although a 20% increase in the desired GDPpc may seem unrealistic, this parameter is not the actual GDPpc, but rather is the expected per capita economic output, which is constrained by resource availability and climate change impacts. On the other hand, for the sensitivity exploration the ranges given to the input parameters are less relevant, since all output perturbations are normalised by the percentage of change of the input parameter with respect to its nominal value.

Uncertainty analysis
The analysis of the 1000 Monte Carlo simulations shows that uncertainty increases over time and the con dence bounds become wider for all output variables analysed. In other words, uncertainty widens as the projected output includes less of the known past. This behaviour is to be partly expectable because some of the tweaked input parameters are constant rates (e.g. growth rate of installed capacity of the different RE technologies).
The EROI of the system disperses from the initial value of ca. 11:1 in 2015 to values between 7:1 to 13.5:1 in 2050 (Fig. 2). Since the EROI of renewable energy generation technologies are nowadays smaller (on average) than those of fossil fuels [38], the lower values observed may correspond to scenarios with higher penetration of RE, therefore scenarios with higher implementation rates of the different RE technologies included in the model. Although the uncertainty for the EROI may seem relatively small, the difference between the highest and the lowest values of the distribution mark the difference between a prosperous society and one at the brink of collapse. Indeed, the lowest EROI values obtained are just above the critical value of 5:1, considered to be the lowest possible value to sustain human societies [43].

Input parameter
these outputs is very large, most of the simulations (con dence bound of 50%) did not reach exponential trends. The GDPpc had more uncertainty than the GDP, because the uncertainty on the population adds to the uncertainty of the GDP itself. In the MEDEAS models, the causation between GDP, energy consumption and CO 2 emissions is bidirectional: more economic activity (measured in GDP) results in more energy consumption, hence higher CO 2 emissions; high atmospheric CO 2 concentrations increase the energy demand for adaptation to climate change and if the energy supply cannot match the increased demand, the GDP is negatively affected. Fig. 3 also shows that none of the biophysical constraints implemented in the model impede the exponential growth of the economy under certain combinations of input parameters. Higher population growth rates are also responsible for higher economic and energetic demand and CO 2 emissions. The uncertainty of the share of electricity generation from RES is very large (Fig. 4), and depends not only on the deployment rates but also on the techno ecological potentials of each RES technology. Also, in a scenario of high population and high economic growth (high energy demand) with less RES penetration (little electricity produced from RES) will have very low shares.
Finally, the large uncertainty of the total land requirements (Fig. 5) comes also from the uncertainty in the penetration rates of each technology, their share (some technologies take more land than others), their techno-ecological potentials and the land available (once land is fully occupied, no more RE infrastructure can be built).

Sensitivity exploration
Very similar results were obtained by using the three different sensitivity measures: the RMSD, the Euclidean distances and the Spider plots. The later measures only a punctual distance at year 2050, but clearly showed the non-linear dependency between inputs and outputs [45] (Fig. 7). The RMSD measures distances for each individual input-output relationship across all simulated times (1995-2050) while the Euclidean distance is the aggregated distance of the vector of all outputs of each simulation, aggregated over time.
Combining the results obtained with the three methodologies, the   most impacting parameters on the model outputs were min_e-nergy_intensity_vs_initial, min_energy_intensity_vs_initial_h, a1_conf_th and variation_nonxdashxenergy_use. Parameters min_energy_intensity_vs_initial and min_e-nergy_intensity_vs_initial_h affect the values of the energy intensities of all economic sectors and households, respectively. The higher their value, the lower potential ef ciency gains, which result in higher TFEC and CO 2 emissions. In addition, the lower the ef ciency, the less energy can be destined to build new RE infrastructure, which results in a reduction of the share of RE electricity. On the other hand, the large impact of a1_conf_th shows the importance of the household transport sector in the whole energy system. To this regard [46] found that a 100% renewable road transport providing the same service as in 2014 would demand 69% less energy.
Finally, the actual demand for fossil fuels for non-energy purposes, described in the model with variable variation_non-energy_use, is poorly understood [47]. Following [47], the demand for each nal fuel (liquids, gases and solids) is endogenously calculated as a linear function of the historic GDP, although this assumption is very uncertain. Therefore, the high impact of this variable on the outputs comes as a result of its structural uncertainty, rather than parametric uncertainty [20].
Based on these results, these 4 parameters will be key when trying to t simulation results between the MEDEAS and other models. On the contrary, seven of the input parameters were seen to have little impact on the results and will be xated in future GSA on the MEDEAS models.
The OAT analysis, although still one of the most widely used techniques for sensitivity analysis in many domains due to its simplicity and low computational cost, is known to have drawbacks when applied to non-linear models [27]. In addition, when using OAT technique a large fraction of the input parameter space is left unanalysed and the effects of input parameters interactions cannot be evaluated [22,48]. Acknowledging these weaknesses, in this work this methodology was used as a previous step, to discard those parameters with negligible effects on the model outputs in a subsequent global sensitivity analysis.

Conclusions
Uncertainty is at the origin of the main criticisms made to IAM models, and unless adequately assessed, it may hold back their use for policy design. Thus, uncertainty (UA) and sensitivity analysis (SA) are key to give robustness to IAMs outputs.
The current work complements the previous work done in Ref. [37], in which the pymedeas models were introduced and an exploratory analysis was performed to assess the impacts of changing model hypotheses. The UA performed here assesses the robustness of the scenarios chosen in the previous work, while the SA serves to identify the most relevant parameters within each scenario. This methodology (hypotheses-scenarios-UA-SA) gives a useful procedure for future energy system and IAMs models analysis.
Results from the Monte Carlo simulations indicate that uncertainty in the outputs is large and increases over time, especially for scenarios with high economic growth expectations. These results also highlight the dangers on the climate and the environment of sustaining a growing economy with fossil fuels. On the other hand, parameters directly affecting sectorial energy ef ciencies and households' transport energy use had the larger impact on the model outputs. These are known to be key elements to achieve the energy transition, and the results of this work support all efforts being made by European policy-makers to this regard. The use of fossil fuels for non-energetic purposes had also a big impact on model projections, but in this case it was attributed to a structural uncertainty rather than a parametric uncertainty.
This work demonstrates how the uncertainty propagates in highly non-linear models with different feedbacks at different sub-models levels [42]. The analysis performed here also highlights the limitations of the system dynamics approach, which works very well for capturing the inner system non-linear feedbacks, but when dealing with rapidly increasing non-linear interactions, uncertainty becomes a 'demon in the machine' that affects the reliability of the model projections [49]. In this sense this work contributes to the analysis of limitations of complex system dynamics models and IAMs [50] in general and in particular to the MEDEAS models.
Finally, the complex interactions between different parameters and strong nonlinearities between inputs and outputs found in this work justify a subsequent and more computationally intensive GSA, that will bene t from the factor prioritization made here.

Declaration of competing interest
The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.