The analyzed product system is constituted by a fleet of BEV passenger cars, which are operated in Italy across their useful life. Plug-in Hybrid Electric Vehicles (PHEV) and hybrid electric vehicles are thus excluded from the study. The functional unit of the study is related to the marginal provision of 1 km transport service by BEV passenger cars. As a starting point, the generic cradle-to-grave dataset “transport, passenger car, electric” from the ecoinvent 3.6 consequential database was considered. Within such dataset, the process related to electricity consumed in BEV use phase constitutes the foreground system of this study, which is refined according to the ESM results. To enhance the interpretation of results, a specific focus is devoted to an intermediate functional unit, as made in Ricardo (2020): i.e., the electricity consumed in the use phase of such system. Firstly, our study analyzes the provision of 1 kWh of low voltage electricity in Italy. Secondly, the supply of the yearly electricity demand associated to the entire BEVs fleet is discussed. Thirdly, the variations of the distributions of hourly LCA impacts associated to the provision of 1 kWh of electricity are shown. Eventually, the variations of LCA results of 1 km transport by BEV, due to the explicit foreground LCI model of use-phase electricity are discussed.

The ESM is chosen coherently with the aims of the study: EnergyPLAN employs an hourly time resolution over a year, while keeping a relatively simple and aggregated structure. The model is able to integrate projections and historical data collected and developed by national research institutes, such as Terna and GSE in the Italian case. Moreover across this study, indications on the Italian national energy strategy are drawn from the “Integrated National Energy and Climate Plan,” which is the latest set of policies submitted to the European Commission in 2019 (Piano Nazionale Integrato Energia e Clima, PNIEC) (MiSE et sl., 2018). The ESM analyses are developed starting from Marmiroli's work (2020).

Concerning LCI modeling, the consequential, long-term approach is adopted. In the foreground system, the electricity process is analyzed in terms of its structural response to changes, being it modeled as a mix of long-term marginal processes. According to the goal of the study, a BAU scenario for 2030 in Italy must be chosen. The most recent studies at the Italian level outline that the current market trend is projected to lead to a range of 1.8–2 million EVs (Enel European House Ambrosetti, 2017; Energy Strategy group, 2018). At the policy level, the PNIEC considers incentives for spreading low-emission passenger vehicles, yet without providing any specific target on expected shares of passenger EVs. Since our analysis only includes BEVs, which do not constitute the totality of EVs, the lowest value of 2030 EVs penetration rate, according to the presented projections, is chosen to represent the BAU scenario. Therefore, we assume that PNIEC objectives include the diffusion of 1.8 million BEVs, which constitutes the BAU scenario. Data on electricity production by source for 2016 and 2030 BAU scenarios are, respectively, collected in accordance to the PNIEC.

Concerning product system boundaries, the useful life of the car modeled in the cited ecoinvent dataset is deducted to be 150,000 km, according to the indications on car weight. This dataset has a global geographical coverage, meaning that its LCI processes are determined through the use of global market mixes. In the present analysis, the geographical coverage of the only use phase is modified, contextualizing it within the Italian country. Within the use phase, the only electricity provided for operating BEVs is included in the foreground system, as it is the only process with both a reasonable dependence on the country in which the vehicle is operated and a noticeable relevance in determining LCA impacts for BEVs. The background system includes production, maintenance, road and wear-related emissions, and End-of-Life (EoL) processes related to BEV; therefore, such system is not changed in this study. Concerning time boundaries, the ecoinvent 3.6 consequential database computes yearly marginal electricity mixes using 2030 as the time horizon year and 2015 as reference year (Vandepaer et al., 2019). The same approach is followed here, yet 2016 is chosen as reference year due to data availability.

Within this framework, five scenarios with different 2030 BEVs diffusion levels into the Italian passenger car fleet are implemented. In accordance to the goal of the study, to show the potential occurrence of non-linearities in LCA results, the penetration of BEVs into the Italian passenger car fleet is modeled up to 100% with a stress test, with 4 steps of 25% each. According to Capros et al. (2016), the size of the European passenger car fleet is not expected to grow from 2016 on. Thus, the absolute level of the Italian passenger car fleet is kept equal to the 2016 value, i.e., 38 million cars (ACI, 2017). Therefore, BEVs constitute 4.75% of the entire fleet in the BAU scenario, while the additional four scenarios model a fleet which is constituted by 25, 50, 75, and 100% BEVs. The last 3 cases are purely speculative, as such values are not reached even in the most optimistic scenarios developed within (Enel European House Ambrosetti, 2017; Energy Strategy group, 2018).

The selected Life Cycle Impact Assessment (LCIA) method is ILCD 2.0 midpoint, as such method is designed for policy decision support (Sala et al., 2012). The selection of the most relevant impact categories for this study is based on three criteria. Firstly, since the focus of the study is related to the Italian marginal electricity mix, the relevance of LCIA results for this process, relative to the cradle-to grave LCIA results, is evaluated. Secondly, the relevance of BEV LCIA results with respect to global impacts is accounted for. Thirdly, a selection within LCIA categories which are similarly influenced by variations of the Italian marginal electricity mix is carried out. Figure 1 focuses on the first two criteria: it shows the normalized BEV LCIA results, grouped by relevant providers, for the “transport, passenger car, electric” ecoinvent dataset, modified using Italian electricity provider for BEV operation. Figure 1 shows that the most relevant LCIA categories for the present analysis can be Non-carcinogenic effects, Carcinogenic effects, Minerals and metals, Climate change total, Freshwater eutrophication, Freshwater ecotoxicity, and Fossils. However, it was observed that LCIA trends for the Fossils and Freshwater ecotoxicity categories are, respectively, similar to Climate change total and Freshwater eutrophication (see Supplementary Figures 7–9), thus according to the third criterion, these categories are not considered. Finally, Respiratory effects, inorganics is included in the study, as it is indicated to be relevant in the Product Environmental Footprint Category Rules (PEFCR) for batteries (Siret et al., 2018), thus it may be important as well for the present study. Selection of impact categories would not have changed if the final results for the electricity dataset were considered: section Additional Discussion of the Supplementary Material provides additional discussion on the selection of impact categories. Moreover, section Outputs of the Supplementary Material provides LCIA results for all the ILCD impact categories.

The present analysis overall requires six ESM simulations to be built: the first is related to historical 2016 data, while the remaining are related to the five 2030 scenarios of BEVs diffusion. Then, marginal electricity mixes for each of the five 2030 scenarios are computed.

Recalling the first aim of the study, three different approaches for computing marginals are here presented, representing three different levels of detail that may be used in LCA for BEV. The first approach computes marginal electricity mixes at the yearly level, while the second one additionally involves the hourly time resolution employed by EnergyPLAN. Finally, the third approach even includes BEV hourly charging profiles. The three approaches are named according to the temporal resolution of marginal electricity productions used to compute them and to the weighting of mixes to account for BEV charging curves, whenever this is performed. Marginal mixes computed through equation 1 are referred to as Yearly. Marginal mixes computed through equation 2 are called Hourly. Eventually, marginal mixes computed through equation 4 are named Hourly weighted, since they account for a weighting process based on BEV charging profiles.

Recalling the second aim of the study, results from the above presented approaches are compared with the conventional linear LCA approach, which is here represented by the use of the ecoinvent 3.6 consequential dataset for the Italian electricity mix. Accordingly, this approach is named Ecoinvent. The ecoinvent long-term consequential system model aims at reflecting long-term consequences of small-scale decisions. Such decisions do not influence key market parameters, thus they are assumed to linearly scale with the size of the related change (Vandepaer et al., 2019). Instead, our analysis investigates the variation of such key market parameters and the consequent occurrence of non-linearities, which are thus quantified in contrast with the Ecoinvent approach. In addition, a comparison with the Ecoinvent approach can outline the relevance of underlying data sources used to compute marginal electricity mixes, as the Italian ecoinvent dataset is derived from European reference scenarios (Vandepaer et al., 2019), differently from our study, as section Modeling of Energy Scenarios outlines.

Eventually, given the third aim of the study, the differences between the 3 BEV-dependent marginal mixes are outlined, too.

The first approach, Yearly, employs Equation 1, presented in Vandepaer et al. (2019), for computing marginal electricity supply shares with respect to 2016, across the five 2030 scenarios:

With i being each of the n unconstrained electricity production technologies, P being the quantity of electricity produced within the related year by technology i, S being the share of technology i in the considered marginal supply mix. Unconstrained electricity production technologies do not include electricity produced in Combined Heat and Power (CHP) mode, as in that case electricity is produced as a secondary product. Whenever the numerator of the equation is negative, S is given the value of 0 (Vandepaer et al., 2019).

The Yearly approach does not consider that even if the eventual P_{i, 2030} is lower than P_{i, 2016}, some technologies may increase their hourly electricity production in 2030, with respect to 2016, across some hours of the year. Thus, stating that a specific electricity production technology is not included in the yearly marginal electricity mix may not be fully true. Section Modeling of Energy Scenarios outlines that for Italy, the 2030 electricity demand is projected to be higher than the 2016 one, but P_{natural gas, 2016} is equal to 126 TWh (GSE, 2018), while the projected P_{natural gas, 2030} is equal to 118 TWh (MiSE MATTM, 2017). Consequently, S_{natural gas, 2030} is null. However, whenever the 2030–2016 increase in electricity demand is not covered by Variable Renewable Electricity Sources (VRES) (i.e., during hours with a low solar irradiation or wind), thermal power plants must cover such increase, thus leading to a positive marginal electricity production within specific hours of the year, as will be shown by Figure 2. The second approach, Hourly, aims at filling this gap using Equation 2:

With h being the considered hour of the year. The same approach was followed by Lund et al. (2010), who computed a long-term yearly average marginal technology based on hourly dynamics. EnergyPLAN employs 8,784 h, even if 2030 is not a leap year. However, whenever hourly distributions are required (e.g., electricity demand, VRES productions), the 2016 distributions are employed for 2030 as well, as a matter of consistency between scenarios. Section Electricity System Specifications and Additional Discussion of the Supplementary Material, respectively, provide explanations and motivations on the use of data distributions within EnergyPLAN.

Still, the Hourly shares do not consider the charging patterns of BEVs, i.e., what the actual electricity consumed by BEVs is. Thus, we propose a third method for calculating marginals, based on the indications of Girardi et al. (2015). Therefore, we introduce W_h, which is defined by Equation 3 and is used to weigh the hourly electricity marginal productions on the charging profiles of BEVs:

With C_h being the amount of electricity which is used to charge the BEVs fleet. The denominator of the equation represents the total hourly marginal electricity production. Then, the Hourly, weighted approach employs Equation 4 to obtain aggregated shares which are representative of the actual electricity consumed by our product system:

Finally, in accordance to the ecoinvent database, all the electricity production terms P_i and P_i,h are computed according to the consumption mix of the country. Due to the penetration of VRES, Excess Electricity Production (EEP) may arise in the ESM simulation. EEP must be either exported or curtailed, thus it is not actually consumed within the country. Consequently, it had to be subtracted from the electricity mix. Hence, the theoretical electricity production P_{2030,theoretical} related to each of the VRES technologies is decreased, according to the values of EEP for the related time step (yearly in case of equation 1 and hourly in case of Equations 2 and 4). The subtraction term shown by Equations 5 and 6 is weighted on the relative electricity production of each of these sources in the considered time interval. Electricity production from Concentrated Solar Power (CSP) is decreased, too. This technology has the possibility to store thermal energy and it may theoretically allow for a higher degree of flexibility of electricity production. However, CSP is still not mature in the Italian market, thus as a conservative hypothesis, this technology is assigned no stabilization property of the grid. Equation 5 is performed at the yearly level and applies to the P_i term of equation 1, for both 2016 and 2030 years:

With i ranging from 1 to 5, which is the number of VRES technologies included within such assessment, i.e., onshore and offshore wind, solar photovoltaic (PV) and CSP, run-of-river hydro. Equation 6 is instead performed at the hourly level, according to the EnergyPLAN outputs:

In order to obtain the LCI and LCIA results related to every scenario, the computed long-term electricity production mixes are inputted to the high and low voltage electricity mix datasets for Italy of the ecoinvent 3.6 consequential database, thus changing the related LCI through the brightway2 software (Mutel, 2017).

Within this subsection, the only data that are key to understand the Results and Discussion sections are presented. A full assessment of the entire energy system of Italy is not related to the aims of the paper, since focus is devoted to the only electricity system. On the demand side, we updated the values of total electricity demand and import/exports, together with explicitly modeling BEVs demand. On the supply side, the electricity production data by source are computed. The evolution of heating and cooling supply by means of technologies which run on electricity, e.g., electric heating and cooling or CHP, is not explicitly modeled here: the related data are kept fixed to the 2016 levels and are discussed in section EnergyPlan of the Supplementary Material, together with additional information on EnergyPLAN technical parameters and equations.

Our analysis overall requires six ESM simulations. The outputs of 2016 simulation are used both as reference for calculating marginals and for calibration of EnergyPLAN outputs. Concerning the 2030 scenarios, the installed capacities of different technologies are computed on the BAU projections of electricity production, according to the PNIEC objectives. Such capacities are not changed across 2030 scenarios. This means that reduced RES curtailment and thermal power plants are assumed to cover the additional BEVs electricity demand across scenarios. Section Discussion of Model Assumptions discusses the motivation behind this assumption and related implications.

Additional information used to model the 2030 BAU scenario is based on the projections elaborated by RSE (2017) for the National Energy Strategy (MiSE MATTM, 2017), which was elaborated in 2017 to create a favorable context for the PNIEC adoption. Quantitative targets of PNIEC that are of interest for this study are the phase out of coal by 2025, together with the 2030 share of RES in electricity production and in the supply of transport demand, which must, respectively, reach 55.4 and 21.6%. The PNIEC policies are estimated to lead to a 36% reduction of greenhouse gases emission from car transport, with respect to 2005, which is mostly due to the electrification of this sector, rather than to biofuels penetration (MiSE et sl., 2018).

EnergyPLAN offers the possibility to perform either a technical analysis, which implements different technical simulation strategies, or a market exchange analysis, which works on the optimization of business-economic profits for each plant (Lund and Thellufsen, 2020). Due to the difficulty in predicting energy prices, a technical simulation strategy is set up. A discussion on the robustness of this analysis is provided in section Discussion of Model Assumptions. Required EnergyPLAN inputs are related to energy demands, production capacities, efficiencies and energy sources. The outputs consist of annual energy production balances, which are then used to compute LCI marginals as described in section LCI Methodological Setup.

All electricity data are here expressed in gross terms, since this is the actual electricity produced. Grid electricity losses are accounted for as in the original ecoinvent dataset. Losses and modeling of hydroelectric pump operation are explained in section Inputs of the Supplementary Material.

Figure 2 shows the marginal electricity mixes of the BAU scenario, resulting from the four approaches outlined in section LCI Methodological Setup. The Yearly approach leads to a mix with a strongly lower share of fossil fuels, with respect to the Ecoinvent one. Indeed, neither natural gas nor biomass condensing power plants increase their yearly production with respect to 2016, as can be clearly seen by Table 3. The Yearly mix is therefore dominated by PV and wind, which account for more than 80% of it. Instead, the Hourly mix appears to be more similar to the Ecoinvent one, with a natural gas fueled plants share of 9.50%, which appears because they operate during hours in which RES electricity production cannot satisfy the additional hourly demand. Finally, the Hourly, weighted mix still shows the presence of natural gas, yet with a smaller contribution with respect to the Hourly marginals. The share of PV is instead higher, since BEVs charging is privileged across the central hours of the day, as Figure 4 will further discuss. Furthermore, the share of hydro, dammed is nearly halved with respect to the Hourly approach. Indeed, as section Modeling of Energy Scenarios explained, when EEP arises, EnergyPLAN decreases dammed hydroelectric power production and instead, water storages are charged by means of pumps.

Then, Figure 3 shows the variability of marginal electricity mixes across the analyzed scenarios, for the 3 non-linear approaches followed in this study. The Ecoinvent shares are not shown here, as they do not vary across scenarios. Due to our hypothesis of keeping VRES capacities constant across scenarios, the general trend is related to the increase of natural gas share to match additional BEV electricity demand. Such increase occurs more steeply across the Yearly marginals, in which the share of natural gas increases from 0% in the 25% BEV scenario to 32.79% in the 100% BEV scenario. Indeed, positive EEP amounts of 12.68 and 3.69 TWh appear, respectively, in the BAU and 25% BEV scenarios. Thus, according to our hypotheses, within the BAU and the 50% BEV scenarios, the increase in BEVs electricity demand is partly covered by EEP exploitation and partly by thermal power plants. This is why the Yearly marginal electricity production from natural gas condensing power plants in the 25% BEV scenario is still not positive (see Supplementary Table 4). From the 50% BEV scenario on, the increase in BEVs electricity demand is only covered by thermal power plants (mostly natural gas, with small contributions of biomass, as can be inferred by Table 3 shares) in condensing mode. This results in an increase in the derivative of the curve related to the natural gas share, which occurs from 25% BEV to 75% BEV scenarios. This variation of derivative is much more evident for the Yearly marginals, since the Hourly and Hourly, weighted approaches account for the dynamics of electricity production within a year, which smooth the increase of natural gas electricity production. From 75% BEV to 100% BEV scenarios, the derivative of the curve related to natural gas share decreases again. The Hourly, weighted marginals show the peculiarity of a steep increase of natural gas share already in the 25% BEV scenario and the related derivative does not decrease when moving toward the 100% BEV scenario, in which the highest natural gas share across the analyzed approaches is reached, being it equal to 38.94%.

To facilitate the interpretation of the Hourly, weighted LCI results, Figure 4 outlines the interdependence between the distribution of BEVs charge profiles, which are related to the weights W_h, and the variation of average hourly electricity mixes across scenarios. Supplementary Figures 2, 3 provide the full variation of both quantities across scenarios. Hourly charge values are shown as a fraction of the electricity charged in the related day, to clearly outline the contribution of each hour.

The distribution of 366 hourly charge profiles results from both the transport demand curve, taken from the EnergyPLAN models developed in Connolly et al. (2015) and the choice of Smart Charge simulation, which stimulates charging when EEP is produced. Even if the distribution of transport demand is fixed across scenarios, distributions of hourly electricity consumed for BEVs charging do change, as the absolute value of transport demand changes with the size of the BEVs fleet, too. Indeed, the BAU and 100% BEV scenarios show strongly different distributions. The distribution of BEVs charge profiles in the BAU scenario shows many outliers, which are due to the EnergyPLAN model seeking to fully exploit electricity produced by PV and CSP. Solar energy is highly relevant in the 2030 Italian mix and is more available within the central hours of the day, differently from other RES. EnergyPLAN applies the same logic to wind energy exploitation, which is responsible for outliers that are not within the central hours of the day. In the BAU scenario, BEVs transport demand is very low and EEP is high, therefore charging occurs mainly to reduce EEP. Conversely, moving toward the 100% BEV scenario, the relevance of BEVs on the total electricity demand increases, thus BEVs are charged also within other hours of the day, which are not necessarily those in which EEP may occur. Hence, additional BEVs make the charging curve flatter and looking more similar to the transport demand curve (see Supplementary Figure 2), with all the 100% BEV ordinates being lower than 25% in Figure 4. Moreover, from the 50% BEV scenario on, EEP is null, thus additional electricity demand is mainly covered by natural gas, which starts penetrating even in the central hours of the day. In particular, Figure 4 shows that the penetration of natural gas especially affects early morning and evening hours. Indeed, night shares do not significantly change across scenarios, as the charging process mainly happens within the day. This means that the increase in electricity demand happens within hours in which solar VRES cannot adequately respond. Therefore, the higher the BEVs diffusion, the higher the contribution of thermal power plants in charging additional BEVs. This trend is fully captured by the Hourly, weighted approach only, as confirmed by the related increase in natural gas share across scenarios in Figure 3. Hence, our hypothesis of keeping the same VRES capacities across different scenarios has particularly significant implications especially on the actual electricity consumed by BEVs.

All the trends outlined in section Technology Shares are reflected onto the related LCIA impacts. The unitary impacts per kWh of marginal electricity, related to the BAU scenario mixes outlined by Figure 2, are shown in Table 4. Then, Figure 5 shows the trends of unitary impacts, relative to the values of BAU scenario with the Hourly, weighted approach. Such representation allows to outline which impact categories are mostly affected by the variation of shares across scenarios and methodological approaches. The Ecoinvent trend is constant, as this approach cannot account for interactions between the electricity LCI and the BEVs fleet.

Carcinogenic effects is the least affected impact indicator, with the curves related to all the 4 approaches being confined within the 80–100% interval. The opposite happens for Climate change total indicator, in which the 3 non-linear curves span across an interval of 60–220%. In the BAU scenario, results for the Ecoinvent approach differ by +39.74% for the Climate change total category and by −49.27% for Minerals and metals category, with respect to the BAU, Hourly, weighted value. Thus, employing different approaches for marginal calculations is a significant source of variability of LCIA indicators. These two categories also show the greatest variations of non-linear curves, as they are strongly dependent by the balance between RES and fossil fuels shares within the electricity mix, due to the unitary impact of each technology (see Supplementary Figure 5). Apart from the Climate change total charts, the difference within the 3 non-linear curves decreases when moving toward the 100% BEV scenario. Indeed as Figure 3 showed, the higher the BEVs demand, the higher the penetration of natural gas. Therefore, within the same scenario and across the 3 non-linear approaches, the differences of shares related to other technologies decrease. Apart from Climate change total, impact indicators are dominated by technologies other than natural gas, thus we see that the 3 non-linear curves progressively approach each other. For instance, the Minerals and metals impact indicator is dominated by PV penetration: the lower the differences in PV shares, the lower the differences in Minerals and metals impacts across the 3 non-linear approaches. Indeed, starting from a 50% wide range, the different curves converge into a very narrow range in the 100% BEV scenario.

The Climate change total indicator computed with the Yearly approach is highly sensitive to the increase of BEVs, as a consequence of the highly variable shares outlined by Figure 3. On the other hand, the Hourly curve shows instead the smoothest trends across LCIA categories, while the opposite holds true for the Hourly, weighted results. Regarding this category, the Hourly and Hourly, weighted curves start from very close values in the BAU scenario, despite the latter having a much smaller share of natural gas in the electricity mix, as shown by Figure 2. This happens because the BAU mix for the Hourly, weighted curve also shows a much smaller share of electricity from hydro, dammed with respect to the same scenario for the Hourly curve. Hydro, dammed electricity is one of the cleanest sources of electricity production within this category, with unitary impacts being one order of magnitude lower than electricity from PV (see Supplementary Figure 5). Moreover, the trend of non-linear curves within this impact category reflects the penetration of natural gas within the marginal electricity supply mix, with a steep growth from 25% BEV to 75% BEV scenario, followed by a decreasing derivative within the 75% BEV and 100% BEV scenarios.

It must be remarked that according to our model, Climate change total is the only category which is negatively affected by higher penetration of BEVs. However, the relative variations of other categories across the 5 scenarios are much smaller. Within the Climate change category, moving from the BAU to the 100% BEV scenario, all the 3 non-linear curves double their unitary impacts.

Across the 6 impact categories, the ranges defined by the 4 curves in the BAU scenario are strongly different with respect to those related to the 100% BEV scenario. For instance, in the Climate change total category, the BAU scenario values range from 60 to 140%, while the 100% BEV scenario values span across 140–220%. Figure 6 shows how this aspect is translated into the LCIA indicators associated to the total marginal electricity demand evolution across additional fleets of BEVs. Such curves are obtained by multiplying the unitary impacts by the yearly electricity demand per vehicle (see section Demand) and by the number of additional BEVs associated to each scenario. From the comparison with the Ecoinvent values, it is evident how linearly upscaling marginal impacts with the number of circulating BEVs (i.e., what would be done with a conventional LCA employing the Ecoinvent approach, or in general employing the only BAU unitary impact indicators) would lead to significant non-linearities, as unitary impacts are found to vary across scenarios, especially for the Climate change total category. Indeed, the linear Ecoinvent approach leads to an underestimation of 8.35 Mton CO₂-eq in the 100% BEV scenario for the Climate change total category, with respect to the Hourly, weighted results. Similarly, linearly scaling the Hourly, weighted results keeping the BAU unitary impact would lead to an underestimation of 12.45 Mton CO₂-eq with respect to the case in which the interaction with the electricity grid is accounted for. In accordance with Figure 5, non-linearities are especially significant for the Climate change total and Minerals and metals category, while this is not the case for the remaining impact categories. Total impacts for penetration of additional BEVs are null in the BAU, as no additional BEV is present in this scenario.

The last result of this section highlights again the capabilities provided by an ESM. Impacts computed with the mixes derived by Equations 2 and 4 are related to an aggregation at the yearly level of 8,784 hourly mixes. Shares defined with the Hourly approach can be stated to weight the hourly electricity mixes by the related amount of marginal electricity production. That is, the higher the marginal electricity production related to a specific hour, the higher the influence of such hour in determining S_{i, 2030}. Shares defined by the Hourly, weighted approach are instead weighted on the BEVs electricity charge of each hour. Given the LCA methodology structure, the weighting process may be done either with mixes (at the LCI level) or with impacts (at the LCIA level). Figure 7 stresses out that the Hourly and Hourly, weighted curves are the result of both an underlying hourly impact distribution, which varies across scenarios, and an aggregation process, which varies across the two approaches and clearly constitutes a source of deviations between aggregated LCIA results. If different BEVs charge profiles were assumed, not only the aggregation of the Hourly, weighted approach would change, but also the underlying distribution would be different, especially at higher number of circulating BEVs: see section Contextualization Within the Literature.

Distributions of hourly impact indicators for the BAU scenario appear to be asymmetrical, with the presence of outliers in the upper part of the charts. This is due to the high EEP in the BAU scenario, which makes the productions of VRES decrease (see Equations 5 and 6) and is generally linked to the activation of both natural gas, for grid stabilization requirements (see section 1.1.1 of the Supplementary Material) and the hydro pump technology, which makes unitary impacts further increase (see section Inputs of the Supplementary Material). Accordingly, the higher the number of BEVs, the more the distributions are symmetrical and less outliers are shown. Carcinogenic effects show a large number of outliers, due to high unitary impacts of wind and solar technologies (see Supplementary Figure 5), which are also the most relevant VRES across marginal mixes. Respiratory effects, inorganics shows the presence of outliers in the BAU and 100% BEV scenario. For the former scenario, the reason lies behind high presence of PV in the mix, which shows the second higher unitary impact in this category. For the latter scenario, the reason lies behind the presence of biomass in the mix, which shows the highest unitary impact in this category.

Minerals and metals distributions yet remain strongly asymmetrical. Given that impacts within this category are dominated by PV (see Supplementary Figure 5), the shape of this distribution can be inferred to be related to the presence of sunny days. It is reasonable that 50% of the distribution is instead confined within a small interval (lower than 2e-6 kg Sb-eq), because these values are related to night hours.

This subsection shows how the variability of the use phase electricity LCI makes the cradle-to-grave LCIA results vary. As in the previous section, Table 5 shows the unitary impacts associated to the provision of 1 km of transport service by passenger electric car in the 4 BAU scenarios. Then, Figure 8 shows the trends of unitary impacts, relative to the values of BAU scenario with the Hourly, weighted approach.

In addition to what section Impact Analysis showed, Figure 8 outlines that the variability of BEV LCA results depends on the relative importance of the electricity process within each of the analyzed impact categories, which can be inferred from Figure 1. Accordingly, even if the electricity process shows a high variability within the Minerals and metals category, Figure 8 shows that the variation of impacts for 1 km of electric transport is within an interval of 10% across the analyzed curves. Indeed, the contribution to impacts of car and battery production processes dominate within this category. LCIA indicators related to Carcinogenic effects and Freshwater eutrophication show a very low variability, due to both low variability of use phase electricity-related impacts and a low relevance of this process within BEVs results. In particular, Minerals and metals and Freshwater eutrophication show an average 6.1% decrease of the Hourly, weighted unitary impacts, across the two extreme scenarios. Within Non-carcinogenic effects and Respiratory effects, inorganics, the differences between the 3 curves is lower than 5%, as well as the deviation with respect to the BAU scenario. This is definitely not the case for the Climate change total category, in which a range of 90.31–110.70% is defined for the BAU scenario across the 4 curves. Moreover, for the Hourly, weighted curve, the 100% BEV scenario shows a 32.52% higher impact with respect to the BAU scenario. Concerning the Yearly plot, the variation is even higher, spanning from 90.31 to 125.90% across the analyzed scenarios, i.e., +39.40% total increase. Therefore, within this category, moving toward high numbers of circulating BEVs, the electricity process becomes more and more dominant in determining the impacts of 1 km BEV transport (see Supplementary Figure 6). The Hourly, weighted curve shows the highest variations between the BAU and the 100% BEV scenarios, which are on average equal to 9.28% across LCIA categories.

The present study coupled the indications of Girardi et al. (2015), which accounted for the interdependence between EV and electricity grid, with the evaluation of different scenarios of EVs diffusion as made in Bohnes et al. (2017) and Ricardo (2020). Bohnes et al. (2017) did not account for changes in the electricity grid operation due to BEVs penetration into the market. They assumed that potential increases of electricity demand could be counterbalanced by efficiency gains in EV. Moreover, their study is focused on an urban area, while our study is related to the country level. We did not consider counterbalancing factors, as we aimed to model relevant interactions with the electricity system, which may occur as outlined by other studies (Girardi et al., 2015; Ricardo, 2020). Girardi et al. (2015) considered energy prices dynamics within his model, which is not the case in this study. In Ricardo (2020), an attributional approach is followed, yet some consequential elements are modeled, too. These are related to the additional electricity demand from EV, economies of scale and decarbonization of production processes over time. Table 6 provides a summary of similar studies that can be compared with the present one, i.e., analyses related to future BEVs.

In addition to the presented results, Huo et al. (2015) evaluated average fuel-cycle emissions for 2025 BEVs operation, finding values between 0.060 and 0.170 kg CO₂-eq/km, respectively, for China BTR and U.S. NPCC regions. Our results are in line with Girardi et al. (2015), especially for the case in which BEVs are charged by electricity produced for 50% by PV. Using the Italian 2030 marginal mix, Girardi et al. (2015) found a much higher result compared to our BAU scenario. This is due to the assumption of the 2030 model employed by that study, which considers BEVs electricity demand as an additional demand, which will be covered by thermal power plants only, assuming that RES are already fully used. This is true for Bauer et al. (2015) as well, which further outlines how BEVs results are influenced by the employed electricity mix. The results found by Burchart-Korol et al. (2018) are also in line with this considerations, as the Polish 2030 electricity mix is constituted by 80% fossil fuels. Bohnes et al. (2017) found a relatively high value of BEVs Climate change result, which is due to retaining a share of coal in the marginal electricity mix and to the higher electricity demand of BEVs. In Ricardo (2020), a sensitivity analysis on electricity chain assumptions over different time spans is performed. Such analysis shows the lowest value of Table 6, due to the larger time span (2050), hypothesis on decarbonization of production processes and a strongly clean electricity mix. Similarly, Cox et al. (2020) integrated 2040 electricity scenarios from integrated assessment models into the ecoinvent database, assessing variations of both the use phase electricity mix and of production processes. Moreover, such study also employed a Monte Carlo global sensitivity analysis for all performance parameters of the assessed EVs. The same authors already assessed such scenarios within a previous study (Cox et al., 2018) which leads to similar considerations, concerning Climate change impacts of future BEVs, thus it was not added to Table 6. Both the TECH 1.5 and the ClimPol scenarios model the implementation of European climate targets. The variability of our results is correctly lower compared to Bauer et al. (2015), Huo et al. (2015), Burchart-Korol et al. (2018), Cox et al. (2020), and Ricardo (2020), which assessed highly variable electricity mixes.

In Mathiesen et al. (2009) and Lund et al. (2010), it is demonstrated that long-term marginal electricity supply cannot be constituted by a single technology only and the determination of such mix requires ESM simulations. Indeed, they stated that affected technologies depend both on the geographical location and on the temporal resolution of the study, especially at higher shares of fluctuating renewable energy sources, analyzing the Danish electricity system. Our study validates such considerations for the Italian context too, as it shows that noticeable differences between the Yearly and Hourly approaches, which, respectively, neglect and account for hourly electricity system dynamics. Moreover, we also show the relevance of the different approaches on LCA results for BEV transport services, further validating the need to perform complex analyses for product systems in which energy supply is a relevant process, as stated in Mathiesen et al. (2009). Furthermore, our study highlights that the LCI of electricity supply is strongly sensitive to the amount of marginal electricity demand (see the variations across scenarios of BEV penetration), and to employing the latest policy indications (see the variations between Ecoinvent and Yearly approaches). Thus, we agree with (Mathiesen et al., 2009; Lund et al., 2010) on the need to perform a scenario analysis, in order to avoid that the LCA results become obsolete due to changes in marginal electricity supply mix. Finally, the literature suggested that the actual electricity used to charge BEVs, which depends on the charging profile and on the interaction with the grid, must be accounted for in order to provide accurate LCA results (Nordelöf et al., 2014; Girardi et al., 2015; Rangaraju et al., 2015). In our study, the Hourly, weighted approach shows significant variations of LCIA results across scenarios. In our model, such variations are indeed related to the variability of the actual electricity consumed by BEVs across scenarios, which are not fully captured by other approaches. Therefore, we confirm the literature indications, as we find that accounting for BEVs charging profiles constitutes an additional source of non-linearities in LCIA results. Hence, we estimate that Hourly, weighted is the most accurate approach in accounting for non-linearities, too.

The 2030 BAU scenario is modeled accounting for the system response to current policies, detailed in the PNIEC, in accordance to the consequential, long-term point of view (see section Modeling of Energy Scenarios). However, our LCI model relies on the assumption that no additional RES capacity is installed to cope with the additional BEVs demand across the different scenarios. Thus, no response, in terms of additional installed capacity, of the 2030 system with respect to the 4 additional scenarios with different levels of circulating BEVs is modeled. Hence, the system response to BEVs which are additional with respect to the BAU scenario is a short-term one, as such response is only constituted by the variation of electricity production (Marmiroli et al., 2018). Accordingly, LCI differences across scenarios are mostly due to natural gas penetration, as shown by Figure 3. Another way to build the LCI model could have been to increase RES installed capacity, in order to be able to meet PNIEC objectives given the BEVs transport demand increase. However, even in this case, in the BAU scenario, the differences between the 4 analyzed approaches would emerge again, as they are not affected by this hypothesis. Changing this hypothesis, we would observe a lower penetration of natural gas in the LCI of marginal mixes across scenarios, which would in turn affect especially the LCIA results of the Climate change total indicator. Only for this category, we estimate that the variability of LCIA results across scenarios would probably decrease. The results of other LCIA categories are instead mainly driven by single RES technologies (e.g., PV is the main driver of Minerals and metals results). Thus, they would vary according to which RES capacities increase and consequently, several different trends may appear. However, changing our main hypothesis would lead to additional problems.

Firstly, adding further hypotheses on additional RES installed capacity for each of the analyzed 2030 scenarios, apart from the BAU one. This means deciding on which type of RES is installed and it would lead to lower reproducibility. Moreover, capacity may not be strongly increased for all RES, as for instance hydroelectric capacity is already almost saturated in Italy (MiSE et sl., 2018).

Secondly, a technical simulation strategy is chosen within EnergyPLAN, thus excluding the dynamics of energy prices. Moreover, the energy system is modeled as if it constituted a single market zone. This is not true for Italy, which is constituted by 6 market zones (Gestore Mercati Energetici, 2020). Moreover, EnergyPLAN explicitly models neither frequency regulation nor inertia within the energy system and the BEVs fleet is modeled as if it was a single battery (Lund and Thellufsen, 2020). Thus, with the current EnergyPLAN setup, we are not able to distinguish whether BEVs demand is generated in the South or in the North of Italy. Moreover, determining accurate BEVs demand distributions at the national scale is challenging, since the actual location of every BEV matters, in determining whether there are sufficient charging stations in a specific area. Grid stabilization requirements are only managed through a specific parameter (see section Electricity System Specifications of the Supplementary Material). This means that simulating strongly different systems, with respect to the 2030 reference one, with the current set up of EnergyPLAN would not lead in any case to more accurate results.

Accordingly, particular attention must be devoted to the geographical context of the study. The country electricity mix is one of the most impacting parameters in the sensitivity analysis performed by Ricardo (2020), even if differences across countries will be lower over the years due to the progressive decarbonization of electricity systems. Within our study, BEVs are charged during the day to reduce EEP and to better use the available solar energy, as captured by the BAU results of the Hourly, weighted approach. However, this is a country-specific result, as it depends on the projected strong reliance of Italy on solar energy. If BEVs were operated in France, which mostly runs on nuclear power, daily BEV charging is likely to not lead to a higher exploitation of solar energy. Another example is constituted by Denmark, which is expected to increase its RES share mostly by means of wind power (Bohnes et al., 2017), which may be available during the night too, thus it has different dynamics than solar power. The general remark is that BEV charging should enable a higher RES exploitation, i.e., when RES are available in the specific area of the study, if the objective is to reduce BEVs carbon footprint. Concerning the geographical scale of the analysis, it may make sense to model high diffusions of BEVs at a smaller scale than the national one. BEVs have the potential to decrease urban air pollution (IEA, 2020) and several studies focused on BEVs diffusion at the city scale (Bohnes et al., 2017; Jian et al., 2017).

Another aspect of interest is related to the BEV charge profiles. Figure 4 outlined that charging BEVs during hours in which VRES cannot provide additional electricity, has particularly strong implications on natural gas usage, since thermal power plants would be the main source of electricity which can adequately respond to such demand. Indeed, the average electricity mixes related to early morning and evening hours are the most affected by additional electricity demand. If different BEVs charge profiles were assumed, or the Smart Charge simulation was not chosen, not only the aggregation of the Hourly, weighted approach would change, but also the underlying distribution (see Figure 7) would be different, especially at higher number of circulating BEVs. If for instance, BEVs were charged during the night, two main effects would happen. Firstly, natural gas production is likely to further increase, as the most relevant VRES in Italy is solar, which is not available during night hours: this would e.g., move the Climate change total distribution upwards. This would shift upwards the Yearly, Hourly and Hourly, weighted curves. Secondly, the Hourly, weighted curve for the Climate change total category would be pushed upwards even more, as BEVs are likely to be charged during hours in which only natural gas can respond to higher BEVs electricity demand. The latter effect would not be captured by the Hourly curve. This further confirms the capability of the Hourly, weighted approach to fully capture changes on the actual electricity consumed by BEVs and the related non-linearities.

Moreover, even if different technical parameters for BEVs modeling were assumed, qualitative trends outlined by this study would still appear, since the present analysis covered a sufficiently wide of BEVs demand. The reference year of the analysis would not have significantly changed results, too. See section Additional Discussion of the Supplementary Material for a detailed discussion.

Furthermore, the only non-linearity studied here is related to electricity, but there could be many more, such as modeling economies of scale for passenger car production and disposal or decarbonization of materials. The frameworks developed in Bohnes et al. (2017) and Ricardo (2020) comprehensively outline different processes that may be coupled with this analysis.

Finally, given the goal of the present analysis, we did not focus on evaluating LCA results for the BEVs fleet which is expected to penetrate in the BAU scenario. Indeed, such analysis would not require a consequential approach and an attributional electricity mix for this scenario should be employed instead, as stated in Dandres et al. (2017) for a similar study on data centers.

Concerning the ESM, EnergyPLAN implies some limitations on energy system simulation. First of all, a technical representation of electricity production is not what occurs in reality. The electricity market is indeed based on the cost-merit order, that defines which power plants can produce and for what amount of electricity (Brondi et al., 2019). However, given the goal of the study, introducing market prices is not strictly required. Some reasonable hypotheses are encompassed by the technical simulation, such as the priority of dispatch for RES (biomass excluded). Moreover, the Smart Charge of BEVs aimed at reducing EEP would happen if electricity prices were considered, too. The literature indicates that charging strategies aimed at both reducing energy cost and damping out RES intermittency are being developed (Mehrjerdi and Rakhshani, 2019). In addition, EnergyPLAN allows to introduce accurate data regarding hourly electricity production by source, electricity demand and imports/export distributions. The 2016 model is thus fully calibrated on available historical data. The 2030 BAU model can be still regarded accurate enough, since total electricity productions are calibrated on the model discussed in (RSE, 2017) for the Italian energy strategy. The issue is related to possible future deviations of electricity production distributions by source, since the electricity grid is expected to further develop, allowing more flexibility through different technologies for electricity storage and distributed energy production paradigms (MiSE et sl., 2018). However, still EnergyPLAN suites the goal of developing a framework which may then be coupled with more accurate models.

Table 7 summarizes strengths and weaknesses related to each of the 4 approaches shown in the modeling, in order to guide the practitioner in the inclusion of non-linearities in the LCA study. Such indications are firstly provided in the following table and then discussed afterwards.

Firstly, updating the ecoinvent dataset with the latest policy indications strongly influences the results, even if equation 1 is used in both the Yearly and the Ecoinvent approaches. Then, moving from the Ecoinvent to the Hourly, weighted approach increases the accuracy of results, but decreases their reproducibility. Moreover, the high amount of data required increases the variability of results, too. Integrating an ESM within the LCA requires a high number of assumptions, which may be hard to manage depending on the geographical scale and on the chosen model. The higher the geographical scale, the more complete the assessment is, yet the higher is the number of data and assumptions required, too. However, if the geographical scale is smaller, e.g., city or regional level, such complexity may be managed with accurate data and models and the variability of results may be more easily assessed. Our framework can flexibly be coupled with different geographical scales or ESMs: the methodology would not change.