diff --git a/envs/environment.yaml b/envs/environment.yaml index 264a5dbd..ee1d1605 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -20,7 +20,7 @@ dependencies: - openpyxl!=3.1.1 - pycountry - seaborn -- snakemake-minimal>=8.5,<8.6 +- snakemake-minimal>=8.5 - memory_profiler - yaml - pytables diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 8d41e893..a17e031f 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -248,8 +248,8 @@ rule build_solar_thermal_profiles: output: solar_thermal=resources("solar_thermal_{scope}_elec_s{simpl}_{clusters}.nc"), resources: - mem_mb=20000, - threads: 16 + mem_mb=8000, + threads: 1 log: logs("build_solar_thermal_profiles_{scope}_s{simpl}_{clusters}.log"), benchmark: @@ -816,6 +816,34 @@ def input_profile_offwind(w): } +rule build_egs_potentials: + params: + snapshots=config_provider("snapshots"), + sector=config_provider("sector"), + costs=config_provider("costs"), + input: + egs_cost="data/egs_costs.json", + regions=resources("regions_onshore_elec_s{simpl}_{clusters}.geojson"), + air_temperature=( + resources("temp_air_total_elec_s{simpl}_{clusters}.nc") + if config_provider("sector", "enhanced_geothermal", "var_cf") + else [] + ), + output: + egs_potentials=resources("egs_potentials_s{simpl}_{clusters}.csv"), + egs_overlap=resources("egs_overlap_s{simpl}_{clusters}.csv"), + egs_capacity_factors=resources("egs_capacity_factors_s{simpl}_{clusters}.csv"), + threads: 1 + resources: + mem_mb=2000, + log: + logs("build_egs_potentials_s{simpl}_{clusters}.log"), + conda: + "../envs/environment.yaml" + script: + "../scripts/build_egs_potentials.py" + + rule prepare_sector_network: params: time_resolution=config_provider("clustering", "temporal", "resolution_sector"), @@ -931,6 +959,21 @@ rule prepare_sector_network: if config_provider("sector", "solar_thermal")(w) else [] ), + egs_potentials=lambda w: ( + resources("egs_potentials_s{simpl}_{clusters}.csv") + if config_provider("sector", "enhanced_geothermal", "enable")(w) + else [] + ), + egs_overlap=lambda w: ( + resources("egs_overlap_s{simpl}_{clusters}.csv") + if config_provider("sector", "enhanced_geothermal", "enable")(w) + else [] + ), + egs_capacity_factors=lambda w: ( + resources("egs_capacity_factors_s{simpl}_{clusters}.csv") + if config_provider("sector", "enhanced_geothermal", "enable")(w) + else [] + ), output: RESULTS + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", diff --git a/scripts/build_egs_potentials.py b/scripts/build_egs_potentials.py index 9216dce2..0fc1ff8a 100644 --- a/scripts/build_egs_potentials.py +++ b/scripts/build_egs_potentials.py @@ -30,6 +30,9 @@ from shapely.geometry import Polygon def prepare_egs_data(egs_file): + """ + Processes the original .json file EGS data to a more human-readable format. + """ with open(egs_file) as f: jsondata = json.load(f) @@ -86,6 +89,65 @@ def prepare_egs_data(egs_file): return egs_data +def prepare_capex(prepared_data): + """ + The source paper provides only data for year and regions where LCOE < + 100Euro/MWh. However, this implementations starts with the costs for 2020 + for all regions and then adjusts the costs according to the user's chosen + setting in the config file. + + As such, for regions where cost data is available only from, say, + 2035, we need to reverse-engineer the costs for 2020. This is done + in the following (unfortunately verbose) function. + """ + + default_year = 2020 + + # obtains all available CAPEX data + capex_df = pd.DataFrame(columns=prepared_data.keys()) + + for year in capex_df.columns: + + year_data = prepared_data[year].groupby("geometry").mean().reset_index() + + for g in year_data.geometry: + + if not g in year_data.geometry.tolist(): + # weird but apparently necessary + continue + + capex_df.loc[g, year] = year_data.loc[ + year_data.geometry == g, "CAPEX" + ].values[0] + + capex_df = capex_df.loc[:, default_year:] + + # fill up missing values assuming cost reduction factors similar to existing values + for sooner, later in zip(capex_df.columns[::-1][1:], capex_df.columns[::-1]): + + missing_mask = capex_df[sooner].isna() + cr_factor = ( + capex_df.loc[~missing_mask, later] / capex_df.loc[~missing_mask, sooner] + ) + + capex_df.loc[missing_mask, sooner] = ( + capex_df.loc[missing_mask, later] / cr_factor.mean() + ) + + # harmonice capacity and CAPEX + p_nom_max = prepared_data[2050].groupby("geometry")["PowerSust"].mean() + p_nom_max = p_nom_max.loc[p_nom_max > 0] + + capex_df = capex_df.loc[p_nom_max.index] + + data = ( + pd.concat((capex_df[default_year], p_nom_max), axis=1) + .reset_index() + .rename(columns={2020: "CAPEX"}) + ) + return gpd.GeoDataFrame(data, geometry=data.geometry) + + def get_capacity_factors(network_regions_file, air_temperatures_file): """ Performance of EGS is higher for lower temperatures, due to more efficient @@ -146,28 +208,9 @@ if __name__ == "__main__": egs_config = snakemake.params["sector"]["enhanced_geothermal"] costs_config = snakemake.params["costs"] - sustainability_factor = egs_config["sustainability_factor"] - # the share of heat that is replenished from the earth's core. - # we are not constraining ourselves to the sustainable share, but - # inversely apply it to our underlying data, which refers to the - # sustainable heat. Source: Relative magnitude of sustainable heat vs - # nonsustainable heat in the paper "From hot rock to useful energy..." - egs_data = prepare_egs_data(snakemake.input.egs_cost) + egs_data = prepare_capex(egs_data) - if egs_config["optimism"]: - egs_data = egs_data[(year := costs_config["year"])] - logger.info( - f"EGS optimism! Building EGS potentials with costs estimated for {year}." - ) - - else: - egs_data = egs_data[(default_year := 2020)] - logger.info( - f"No EGS optimism! Building EGS potentials with {default_year} costs." - ) - - egs_data = egs_data.loc[egs_data["PowerSust"] > 0].reset_index(drop=True) egs_regions = egs_data.geometry network_regions = ( @@ -188,7 +231,12 @@ if __name__ == "__main__": overlap_matrix.to_csv(snakemake.output["egs_overlap"]) - # consider not only replenished heat + # the share of heat that is replenished from the earth's core. + # we are not constraining ourselves to the sustainable share, but + # inversely apply it to our underlying data, which refers to the + # sustainable heat. Source: Relative magnitude of sustainable heat vs + # nonsustainable heat in the paper "From hot rock to useful energy..." + sustainability_factor = egs_config["sustainability_factor"] egs_data["p_nom_max"] = egs_data["PowerSust"] / sustainability_factor egs_data[["p_nom_max", "CAPEX"]].to_csv(snakemake.output["egs_potentials"]) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 2d895c65..71445735 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -190,7 +190,7 @@ def define_spatial(nodes, options): # deep geothermal spatial.geothermal_heat = SimpleNamespace() - spatial.geothermal_heat.nodes = ["EU deep geothermal"] + spatial.geothermal_heat.nodes = ["EU enhanced geothermal systems"] spatial.geothermal_heat.locations = ["EU"] return spatial @@ -3599,15 +3599,11 @@ def add_enhanced_geothermal( # under egs optimism, the expected cost reductions also cover costs for ORC # hence, the ORC costs are no longer taken from technology-data - orc_capex = ( - costs.at["organic rankine cycle", "investment"] - if not snakemake.params.sector["enhanced_geothermal_optimism"] - else 0.0 - ) + orc_capex = costs.at["organic rankine cycle", "investment"] # cost for ORC is subtracted, as it is already included in the geothermal cost. # The orc cost are attributed to a separate link representing the ORC. - # also capital_cost conversion Eurofalse/kW -> Euro/MW + # also capital_cost conversion Euro/kW -> Euro/MW egs_potentials["capital_cost"] = ( (egs_annuity + FOM / (1.0 + FOM)) @@ -3619,8 +3615,8 @@ def add_enhanced_geothermal( egs_potentials["capital_cost"] > 0 ).all(), "Error in EGS cost, negative values found." - plant_annuity = calculate_annuity(costs.at["organic rankine cycle", "lifetime"], dr) - plant_capital_cost = (plant_annuity + FOM / (1 + FOM)) * orc_capex * Nyears + orc_annuity = calculate_annuity(costs.at["organic rankine cycle", "lifetime"], dr) + orc_capital_cost = (orc_annuity + FOM / (1 + FOM)) * orc_capex * Nyears efficiency_orc = costs.at["organic rankine cycle", "efficiency"] efficiency_dh = costs.at["geothermal", "efficiency residential heat"] @@ -3647,19 +3643,22 @@ def add_enhanced_geothermal( p_nom_extendable=True, ) - if snakemake.params.sector["enhanced_geothermal_var_cf"]: + if egs_config["var_cf"]: efficiency = pd.read_csv( snakemake.input.egs_capacity_factors, parse_dates=True, index_col=0 ) logger.info("Adding Enhanced Geothermal with time-varying capacity factors.") else: - efficiency = pd.Series(1, overlap.index) + efficiency = 1.0 # if urban central heat exists, adds geothermal as CHP as_chp = "urban central heat" in n.loads.carrier.unique() if as_chp: - logger.info("Adding Enhanced Geothermal as Combined Heat and Power.") + logger.info("Adding EGS as Combined Heat and Power.") + + else: + logger.info("Adding EGS for Electricity Only.") for bus, bus_overlap in overlap.iterrows(): if not bus_overlap.sum(): @@ -3674,34 +3673,36 @@ def add_enhanced_geothermal( bus_egs["p_nom_max"] = bus_egs["p_nom_max"].multiply(bus_overlap) bus_egs = bus_egs.loc[bus_egs.p_nom_max > 0.0] - if egs_config["performant"]: - bus_egs = bus_egs.sort_values(by="capital_cost").iloc[:1] - appendix = pd.Index([""]) - else: - appendix = " " + pd.Index(np.arange(len(bus_egs)).astype(str)) + appendix = " " + pd.Index(np.arange(len(bus_egs)).astype(str)) # add surface bus - n.add( + n.madd( "Bus", - f"{bus} geothermal heat surface", + pd.Index([f"{bus} geothermal heat surface"]), + location=bus, + unit="MWh_th", + carrier="geothermal heat", ) bus_egs.index = np.arange(len(bus_egs)).astype(str) well_name = f"{bus} enhanced geothermal" + appendix - bus_eta = pd.concat( - (efficiency[bus].rename(idx) for idx in well_name), - axis=1, - ) + if egs_config["var_cf"]: + bus_eta = pd.concat( + (efficiency[bus].rename(idx) for idx in well_name), + axis=1, + ) + else: + bus_eta = efficiency p_nom_max = bus_egs["p_nom_max"] capital_cost = bus_egs["capital_cost"] bus1 = pd.Series(f"{bus} geothermal heat surface", well_name) + # adding geothermal wells as multiple generators to represent supply curve n.madd( "Link", well_name, - location=bus, bus0=spatial.geothermal_heat.nodes, bus1=bus1, carrier="geothermal heat", @@ -3711,6 +3712,7 @@ def add_enhanced_geothermal( efficiency=bus_eta, ) + # adding Organic Rankine Cycle as a single link n.add( "Link", bus + " geothermal organic rankine cycle", @@ -3718,7 +3720,7 @@ def add_enhanced_geothermal( bus1=bus, p_nom_extendable=True, carrier="geothermal organic rankine cycle", - capital_cost=plant_capital_cost * efficiency_orc, + capital_cost=orc_capital_cost * efficiency_orc, efficiency=efficiency_orc, ) @@ -3729,7 +3731,7 @@ def add_enhanced_geothermal( bus0=f"{bus} geothermal heat surface", bus1=bus + " urban central heat", carrier="geothermal district heat", - capital_cost=plant_capital_cost + capital_cost=orc_capital_cost * efficiency_orc * costs.at["geothermal", "district heating cost"], efficiency=efficiency_dh, @@ -3745,8 +3747,8 @@ def add_enhanced_geothermal( # Hence, it is counter-intuitive to install it at the surface bus, # this is however the more lean and computationally efficient solution. - max_hours = egs_config["reservoir_max_hours"] - boost = egs_config["reservoir_max_boost"] + max_hours = egs_config["max_hours"] + boost = egs_config["max_boost"] n.add( "StorageUnit", @@ -3756,6 +3758,7 @@ def add_enhanced_geothermal( p_nom_extendable=True, p_min_pu=-boost, max_hours=max_hours, + cyclic_state_of_charge=True, ) @@ -3884,6 +3887,12 @@ if __name__ == "__main__": if options["electricity_distribution_grid"]: insert_electricity_distribution_grid(n, costs) + if options["enhanced_geothermal"].get("enable", False): + logger.info("Adding Enhanced Geothermal Systems (EGS).") + add_enhanced_geothermal( + n, snakemake.input["egs_potentials"], snakemake.input["egs_overlap"], costs + ) + maybe_adjust_costs_and_potentials(n, snakemake.params["adjustments"]) if options["gas_distribution_grid"]: @@ -3911,12 +3920,6 @@ if __name__ == "__main__": if options.get("cluster_heat_buses", False) and not first_year_myopic: cluster_heat_buses(n) - if options["enhanced_geothermal"].get("enable", False): - logger.info("Adding Enhanced Geothermal Potential.") - add_enhanced_geothermal( - n, snakemake.input["egs_potentials"], snakemake.input["egs_overlap"], costs - ) - n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) sanitize_carriers(n, snakemake.config) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 83c23691..9477ce6d 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -801,7 +801,7 @@ def add_geothermal_chp_constraint(n): elec_index = n.links.loc[ n.links.carrier == "geothermal organic rankine cycle" ].index - heat_index = n.links.loc[n.links.carrier == "geothermal heat district heat"].index + heat_index = n.links.loc[n.links.carrier == "geothermal district heat"].index p_nom_lhs = ( n.model["Link-p_nom"].loc[heat_index] - n.model["Link-p_nom"].loc[elec_index] @@ -824,7 +824,7 @@ def add_flexible_egs_constraint(n): n.model.add_constraints( p_nom_lhs <= p_nom_rhs, - name="Upper bounds the charging capacity of the storage unit", + name="Upper bounds the charging capacity of the geothermal reservoir according to the well capacity", ) @@ -954,28 +954,6 @@ def solve_network(n, config, solving, **kwargs): n.model.print_infeasibilities() raise RuntimeError("Solving status 'infeasible'") - # check if enhanced_geothermal_performant might have changed model results - if ( - snakemake.config["sector"]["enhanced_geothermal"] - and snakemake.config["sector"]["enhanced_geothermal_performant"] - ): - mask = (mask := n.links.carrier == "geothermal heat") & ( - n.links.loc[mask, "p_nom_max"] > 0.0 - ) - - saturated = n.links.loc[mask, "p_nom_max"] == n.links.loc[mask, "p_nom_opt"] - - if saturated.any(): - logger.warning( - ( - "Potential for enhanced geothermal heat is saturated at bus(es):\n" - f"{', '.join(n.links.loc[saturated.loc[saturated.astype(bool)].index, 'location'].tolist())}.\n" - "Consider setting config['sector']['enhanced_geothermal_performant'] to False." - ) - ) - - return n - if __name__ == "__main__": if "snakemake" not in globals():