diff --git a/config/config.default.yaml b/config/config.default.yaml
index 7dc0cf76f..6d2ebd9f9 100644
--- a/config/config.default.yaml
+++ b/config/config.default.yaml
@@ -454,6 +454,8 @@ sector:
   hydrogen_turbine: false
   SMR: true
   SMR_cc: true
+  regional_methanol_demand: false  #set to true if regional CO2 constraints needed
+  regional_oil_demand: false  #set to true if regional CO2 constraints needed
   regional_co2_sequestration_potential:
     enable: false
     attribute: 'conservative estimate Mt'
@@ -799,6 +801,7 @@ plotting:
     Coal: '#545454'
     coal: '#545454'
     Coal marginal: '#545454'
+    coal for industry: '#343434'
     solid: '#545454'
     Lignite: '#826837'
     lignite: '#826837'
diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk
index 8a93d24a8..06fd9b79a 100644
--- a/rules/solve_myopic.smk
+++ b/rules/solve_myopic.smk
@@ -88,11 +88,13 @@ rule solve_sector_network_myopic:
         co2_sequestration_potential=config["sector"].get(
             "co2_sequestration_potential", 200
         ),
+        countries=config["countries"],
     input:
         network=RESULTS
         + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
         costs="data/costs_{planning_horizons}.csv",
         config=RESULTS + "config.yaml",
+        co2_totals_name=RESOURCES + "co2_totals.csv",
     output:
         RESULTS
         + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
diff --git a/scripts/add_brownfield.py b/scripts/add_brownfield.py
index 741025801..ffdaf46be 100644
--- a/scripts/add_brownfield.py
+++ b/scripts/add_brownfield.py
@@ -120,6 +120,39 @@ def add_brownfield(n, n_p, year):
             n.links.loc[new_pipes, "p_nom_min"] = 0.0
 
 
+def disable_grid_expansion_if_LV_limit_hit(n):
+    if not "lv_limit" in n.global_constraints.index:
+        return
+
+    # calculate minimum LV
+    attr = "nom_min"
+    dc = n.links.index[n.links.carrier == "DC"]
+    tot = (n.lines["s_" + attr] * n.lines["length"]).sum() + (
+        n.links.loc[dc, "p_" + attr] * n.links.loc[dc, "length"]
+    ).sum()
+
+    diff = n.global_constraints.at["lv_limit", "constant"] - tot
+
+    # allow small numerical differences
+    limit = 1
+
+    if diff < limit:
+        logger.info(
+            f"LV is already reached (gap {diff}), disabling expansion and LV limit"
+        )
+        expandable_acs = n.lines.index[n.lines.s_nom_extendable]
+        n.lines.loc[expandable_acs, "s_nom_extendable"] = False
+        n.lines.loc[expandable_acs, "s_nom"] = n.lines.loc[expandable_acs, "s_nom_min"]
+
+        expandable_dcs = n.links.index[
+            n.links.p_nom_extendable & (n.links.carrier == "DC")
+        ]
+        n.links.loc[expandable_dcs, "p_nom_extendable"] = False
+        n.links.loc[expandable_dcs, "p_nom"] = n.links.loc[expandable_dcs, "p_nom_min"]
+
+        n.global_constraints.drop("lv_limit", inplace=True)
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -150,5 +183,7 @@ def add_brownfield(n, n_p, year):
 
     add_brownfield(n, n_p, year)
 
+    disable_grid_expansion_if_LV_limit_hit(n)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output[0])
diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py
index 1474b0042..7ddc6b1da 100644
--- a/scripts/add_existing_baseyear.py
+++ b/scripts/add_existing_baseyear.py
@@ -303,7 +303,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
         else:
             bus0 = vars(spatial)[carrier[generator]].nodes
             if "EU" not in vars(spatial)[carrier[generator]].locations:
-                bus0 = bus0.intersection(capacity.index + " gas")
+                bus0 = bus0.intersection(capacity.index + " " + carrier[generator])
 
             # check for missing bus
             missing_bus = pd.Index(bus0).difference(n.buses.index)
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index ee2f0e3c6..f1ddce2d6 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -127,15 +127,43 @@ def define_spatial(nodes, options):
     spatial.h2.locations = nodes
 
     # methanol
+
+    # beware: unlike other carriers, uses locations rather than locations+carriername
+    # this allows to avoid separation between nodes and locations
+
     spatial.methanol = SimpleNamespace()
+
     spatial.methanol.nodes = ["EU methanol"]
     spatial.methanol.locations = ["EU"]
 
+    if options["regional_methanol_demand"]:
+        spatial.methanol.demand_locations = nodes
+        spatial.methanol.shipping = nodes + " shipping methanol"
+    else:
+        spatial.methanol.demand_locations = ["EU"]
+        spatial.methanol.shipping = ["EU shipping methanol"]
+
     # oil
     spatial.oil = SimpleNamespace()
+
     spatial.oil.nodes = ["EU oil"]
     spatial.oil.locations = ["EU"]
 
+    if options["regional_oil_demand"]:
+        spatial.oil.demand_locations = nodes
+        spatial.oil.naphtha = nodes + " naphtha for industry"
+        spatial.oil.kerosene = nodes + " kerosene for aviation"
+        spatial.oil.shipping = nodes + " shipping oil"
+        spatial.oil.agriculture_machinery = nodes + " agriculture machinery oil"
+        spatial.oil.land_transport = nodes + " land transport oil"
+    else:
+        spatial.oil.demand_locations = ["EU"]
+        spatial.oil.naphtha = ["EU naphtha for industry"]
+        spatial.oil.kerosene = ["EU kerosene for aviation"]
+        spatial.oil.shipping = ["EU shipping oil"]
+        spatial.oil.agriculture_machinery = ["EU agriculture machinery oil"]
+        spatial.oil.land_transport = ["EU land transport oil"]
+
     # uranium
     spatial.uranium = SimpleNamespace()
     spatial.uranium.nodes = ["EU uranium"]
@@ -1467,8 +1495,8 @@ def add_land_transport(n, costs):
         n.madd(
             "Bus",
             nodes,
-            location=nodes,
             suffix=" EV battery",
+            location=nodes,
             carrier="Li ion",
             unit="MWh_el",
         )
@@ -1567,29 +1595,42 @@ def add_land_transport(n, costs):
 
         ice_efficiency = options["transport_internal_combustion_efficiency"]
 
+        p_set_land_transport_oil = (
+            ice_share
+            / ice_efficiency
+            * transport[nodes].rename(columns=lambda x: x + " land transport oil")
+        )
+
+        if not options["regional_oil_demand"]:
+            p_set_land_transport_oil = p_set_land_transport_oil.sum(axis=1).to_frame(
+                name="EU land transport oil"
+            )
+
         n.madd(
-            "Load",
-            nodes,
-            suffix=" land transport oil",
-            bus=spatial.oil.nodes,
+            "Bus",
+            spatial.oil.land_transport,
+            location=spatial.oil.demand_locations,
             carrier="land transport oil",
-            p_set=ice_share / ice_efficiency * transport[nodes],
+            unit="land transport",
         )
 
-        co2 = (
-            ice_share
-            / ice_efficiency
-            * transport[nodes].sum().sum()
-            / nhours
-            * costs.at["oil", "CO2 intensity"]
+        n.madd(
+            "Load",
+            spatial.oil.land_transport,
+            bus=spatial.oil.land_transport,
+            carrier="land transport oil",
+            p_set=p_set_land_transport_oil,
         )
 
-        n.add(
-            "Load",
-            "land transport oil emissions",
-            bus="co2 atmosphere",
-            carrier="land transport oil emissions",
-            p_set=-co2,
+        n.madd(
+            "Link",
+            spatial.oil.land_transport,
+            bus0=spatial.oil.nodes,
+            bus1=spatial.oil.land_transport,
+            bus2="co2 atmosphere",
+            carrier="land transport oil",
+            efficiency2=costs.at["oil", "CO2 intensity"],
+            p_nom_extendable=True,
         )
 
 
@@ -2419,9 +2460,14 @@ def add_industry(n, costs):
         efficiency=1.0,
     )
 
+    if len(spatial.biomass.industry_cc) <= 1 and len(spatial.co2.nodes) > 1:
+        link_names = nodes + " " + spatial.biomass.industry_cc
+    else:
+        link_names = spatial.biomass.industry_cc
+
     n.madd(
         "Link",
-        spatial.biomass.industry_cc,
+        link_names,
         bus0=spatial.biomass.nodes,
         bus1=spatial.biomass.industry,
         bus2="co2 atmosphere",
@@ -2609,48 +2655,44 @@ def add_industry(n, costs):
         efficiency = (
             options["shipping_oil_efficiency"] / options["shipping_methanol_efficiency"]
         )
-        p_set_methanol = shipping_methanol_share * p_set.sum() * efficiency
 
-        n.madd(
-            "Load",
-            spatial.methanol.nodes,
-            suffix=" shipping methanol",
-            bus=spatial.methanol.nodes,
-            carrier="shipping methanol",
-            p_set=p_set_methanol,
+        p_set_methanol = (
+            shipping_methanol_share
+            * p_set.rename(lambda x: x + " shipping methanol")
+            * efficiency
         )
 
-        # CO2 intensity methanol based on stoichiometric calculation with 22.7 GJ/t methanol (32 g/mol), CO2 (44 g/mol), 277.78 MWh/TJ = 0.218 t/MWh
-        co2 = p_set_methanol / options["MWh_MeOH_per_tCO2"]
+        if not options["regional_methanol_demand"]:
+            p_set_methanol = p_set_methanol.sum()
 
-        n.add(
-            "Load",
-            "shipping methanol emissions",
-            bus="co2 atmosphere",
-            carrier="shipping methanol emissions",
-            p_set=-co2,
+        n.madd(
+            "Bus",
+            spatial.methanol.shipping,
+            location=spatial.methanol.demand_locations,
+            carrier="shipping methanol",
+            unit="MWh_LHV",
         )
 
-    if shipping_oil_share:
-        p_set_oil = shipping_oil_share * p_set.sum()
-
         n.madd(
             "Load",
-            spatial.oil.nodes,
-            suffix=" shipping oil",
-            bus=spatial.oil.nodes,
-            carrier="shipping oil",
-            p_set=p_set_oil,
+            spatial.methanol.shipping,
+            bus=spatial.methanol.shipping,
+            carrier="shipping methanol",
+            p_set=p_set_methanol,
         )
 
-        co2 = p_set_oil * costs.at["oil", "CO2 intensity"]
-
-        n.add(
-            "Load",
-            "shipping oil emissions",
-            bus="co2 atmosphere",
-            carrier="shipping oil emissions",
-            p_set=-co2,
+        n.madd(
+            "Link",
+            spatial.methanol.shipping,
+            bus0=spatial.methanol.nodes,
+            bus1=spatial.methanol.shipping,
+            bus2="co2 atmosphere",
+            carrier="shipping methanol",
+            p_nom_extendable=True,
+            efficiency2=1
+            / options[
+                "MWh_MeOH_per_tCO2"
+            ],  # CO2 intensity methanol based on stoichiometric calculation with 22.7 GJ/t methanol (32 g/mol), CO2 (44 g/mol), 277.78 MWh/TJ = 0.218 t/MWh
         )
 
     if "oil" not in n.buses.carrier.unique():
@@ -2666,7 +2708,8 @@ def add_industry(n, costs):
         # could correct to e.g. 0.001 EUR/kWh * annuity and O&M
         n.madd(
             "Store",
-            [oil_bus + " Store" for oil_bus in spatial.oil.nodes],
+            spatial.oil.nodes,
+            suffix=" Store",
             bus=spatial.oil.nodes,
             e_nom_extendable=True,
             e_cyclic=True,
@@ -2683,6 +2726,39 @@ def add_industry(n, costs):
             marginal_cost=costs.at["oil", "fuel"],
         )
 
+    if shipping_oil_share:
+        p_set_oil = shipping_oil_share * p_set.rename(lambda x: x + " shipping oil")
+
+        if not options["regional_oil_demand"]:
+            p_set_oil = p_set_oil.sum()
+
+        n.madd(
+            "Bus",
+            spatial.oil.shipping,
+            location=spatial.oil.demand_locations,
+            carrier="shipping oil",
+            unit="MWh_LHV",
+        )
+
+        n.madd(
+            "Load",
+            spatial.oil.shipping,
+            bus=spatial.oil.shipping,
+            carrier="shipping oil",
+            p_set=p_set_oil,
+        )
+
+        n.madd(
+            "Link",
+            spatial.oil.shipping,
+            bus0=spatial.oil.nodes,
+            bus1=spatial.oil.shipping,
+            bus2="co2 atmosphere",
+            carrier="shipping oil",
+            p_nom_extendable=True,
+            efficiency2=costs.at["oil", "CO2 intensity"],
+        )
+
     if options["oil_boilers"]:
         nodes_heat = create_nodes_for_heat_sector()[0]
 
@@ -2724,53 +2800,101 @@ def add_industry(n, costs):
         lifetime=costs.at["Fischer-Tropsch", "lifetime"],
     )
 
+    # naphtha
     demand_factor = options.get("HVC_demand_factor", 1)
-    p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nhours
     if demand_factor != 1:
         logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.")
 
+    p_set_plastics = (
+        demand_factor
+        * industrial_demand.loc[nodes, "naphtha"].rename(
+            lambda x: x + " naphtha for industry"
+        )
+        / nhours
+    )
+
+    if not options["regional_oil_demand"]:
+        p_set_plastics = p_set_plastics.sum()
+
+    n.madd(
+        "Bus",
+        spatial.oil.naphtha,
+        location=spatial.oil.demand_locations,
+        carrier="naphtha for industry",
+        unit="MWh_LHV",
+    )
+
     n.madd(
         "Load",
-        ["naphtha for industry"],
-        bus=spatial.oil.nodes,
+        spatial.oil.naphtha,
+        bus=spatial.oil.naphtha,
         carrier="naphtha for industry",
-        p_set=p_set,
+        p_set=p_set_plastics,
     )
 
+    # some CO2 from naphtha are process emissions from steam cracker
+    # rest of CO2 released to atmosphere either in waste-to-energy or decay
+    process_co2_per_naphtha = (
+        industrial_demand.loc[nodes, "process emission from feedstock"].sum()
+        / industrial_demand.loc[nodes, "naphtha"].sum()
+    )
+    emitted_co2_per_naphtha = costs.at["oil", "CO2 intensity"] - process_co2_per_naphtha
+
+    n.madd(
+        "Link",
+        spatial.oil.naphtha,
+        bus0=spatial.oil.nodes,
+        bus1=spatial.oil.naphtha,
+        bus2="co2 atmosphere",
+        bus3=spatial.co2.process_emissions,
+        carrier="naphtha for industry",
+        p_nom_extendable=True,
+        efficiency2=emitted_co2_per_naphtha,
+        efficiency3=process_co2_per_naphtha,
+    )
+
+    # aviation
     demand_factor = options.get("aviation_demand_factor", 1)
+    if demand_factor != 1:
+        logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.")
+
     all_aviation = ["total international aviation", "total domestic aviation"]
+
     p_set = (
         demand_factor
-        * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum()
+        * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1)
         * 1e6
         / nhours
+    ).rename(lambda x: x + " kerosene for aviation")
+
+    if not options["regional_oil_demand"]:
+        p_set = p_set.sum()
+
+    n.madd(
+        "Bus",
+        spatial.oil.kerosene,
+        location=spatial.oil.demand_locations,
+        carrier="kerosene for aviation",
+        unit="MWh_LHV",
     )
-    if demand_factor != 1:
-        logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.")
 
     n.madd(
         "Load",
-        ["kerosene for aviation"],
-        bus=spatial.oil.nodes,
+        spatial.oil.kerosene,
+        bus=spatial.oil.kerosene,
         carrier="kerosene for aviation",
         p_set=p_set,
     )
 
-    # NB: CO2 gets released again to atmosphere when plastics decay or kerosene is burned
-    # except for the process emissions when naphtha is used for petrochemicals, which can be captured with other industry process emissions
-    # tco2 per hour
-    co2_release = ["naphtha for industry", "kerosene for aviation"]
-    co2 = (
-        n.loads.loc[co2_release, "p_set"].sum() * costs.at["oil", "CO2 intensity"]
-        - industrial_demand.loc[nodes, "process emission from feedstock"].sum() / nhours
-    )
-
-    n.add(
-        "Load",
-        "oil emissions",
-        bus="co2 atmosphere",
-        carrier="oil emissions",
-        p_set=-co2,
+    n.madd(
+        "Link",
+        spatial.oil.kerosene,
+        bus0=spatial.oil.nodes,
+        bus1=spatial.oil.kerosene,
+        bus2="co2 atmosphere",
+        carrier="kerosene for aviation",
+        p_nom_extendable=True,
+        efficiency2=costs.at["oil", "CO2 intensity"],
     )
 
     # TODO simplify bus expression
@@ -2821,7 +2945,7 @@ def add_industry(n, costs):
         unit="t_co2",
     )
 
-    sel = ["process emission", "process emission from feedstock"]
+    sel = ["process emission"]
     if options["co2_spatial"] or options["co2network"]:
         p_set = (
             -industrial_demand.loc[nodes, sel]
@@ -2832,8 +2956,6 @@ def add_industry(n, costs):
     else:
         p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / nhours
 
-    # this should be process emissions fossil+feedstock
-    # then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand
     n.madd(
         "Load",
         spatial.co2.process_emissions,
@@ -2994,9 +3116,9 @@ def add_agriculture(n, costs):
             f"Total agriculture machinery shares sum up to {total_share:.2%}, corresponding to increased or decreased demand assumptions."
         )
 
-    machinery_nodal_energy = pop_weighted_energy_totals.loc[
-        nodes, "total agriculture machinery"
-    ]
+    machinery_nodal_energy = (
+        pop_weighted_energy_totals.loc[nodes, "total agriculture machinery"] * 1e6
+    )
 
     if electric_share > 0:
         efficiency_gain = (
@@ -3010,36 +3132,44 @@ def add_agriculture(n, costs):
             suffix=" agriculture machinery electric",
             bus=nodes,
             carrier="agriculture machinery electric",
-            p_set=electric_share
-            / efficiency_gain
-            * machinery_nodal_energy
-            * 1e6
-            / nhours,
+            p_set=electric_share / efficiency_gain * machinery_nodal_energy / nhours,
         )
 
     if oil_share > 0:
+        p_set = (
+            oil_share
+            * machinery_nodal_energy.rename(lambda x: x + " agriculture machinery oil")
+            / nhours
+        )
+
+        if not options["regional_oil_demand"]:
+            p_set = p_set.sum()
+
         n.madd(
-            "Load",
-            ["agriculture machinery oil"],
-            bus=spatial.oil.nodes,
+            "Bus",
+            spatial.oil.agriculture_machinery,
+            location=spatial.oil.demand_locations,
             carrier="agriculture machinery oil",
-            p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nhours,
+            unit="MWh_LHV",
         )
 
-        co2 = (
-            oil_share
-            * machinery_nodal_energy.sum()
-            * 1e6
-            / nhours
-            * costs.at["oil", "CO2 intensity"]
+        n.madd(
+            "Load",
+            spatial.oil.agriculture_machinery,
+            bus=spatial.oil.agriculture_machinery,
+            carrier="agriculture machinery oil",
+            p_set=p_set,
         )
 
-        n.add(
-            "Load",
-            "agriculture machinery oil emissions",
-            bus="co2 atmosphere",
-            carrier="agriculture machinery oil emissions",
-            p_set=-co2,
+        n.madd(
+            "Link",
+            spatial.oil.agriculture_machinery,
+            bus0=spatial.oil.nodes,
+            bus1=spatial.oil.agriculture_machinery,
+            bus2="co2 atmosphere",
+            carrier="agriculture machinery oil",
+            p_nom_extendable=True,
+            efficiency2=costs.at["oil", "CO2 intensity"],
         )
 
 
diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index 224d47147..4bdbb543f 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -350,7 +350,7 @@ def prepare_network(
         # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full
         # TODO: retrieve color and nice name from config
         n.add("Carrier", "load", color="#dd2e23", nice_name="Load shedding")
-        buses_i = n.buses.query("carrier == 'AC'").index
+        buses_i = n.buses.index
         if not np.isscalar(load_shedding):
             # TODO: do not scale via sign attribute (use Eur/MWh instead of Eur/kWh)
             load_shedding = 1e2  # Eur/kWh
@@ -792,6 +792,20 @@ def extra_functionality(n, snapshots):
         add_carbon_budget_constraint(n, snapshots)
         add_retrofit_gas_boiler_constraint(n, snapshots)
 
+    if "additional_functionality" in snakemake.input.keys():
+        import importlib
+        import os
+        import sys
+
+        sys.path.append(os.path.dirname(snakemake.input.additional_functionality))
+        additional_functionality = importlib.import_module(
+            os.path.splitext(
+                os.path.basename(snakemake.input.additional_functionality)
+            )[0]
+        )
+
+        additional_functionality.additional_functionality(n, snapshots, snakemake)
+
 
 def solve_network(n, config, solving, opts="", **kwargs):
     set_of_options = solving["solver"]["options"]
@@ -839,6 +853,10 @@ def solve_network(n, config, solving, opts="", **kwargs):
             f"Solving status '{status}' with termination condition '{condition}'"
         )
     if "infeasible" in condition:
+        m = n.model
+        labels = m.compute_infeasibilities()
+        print(labels)
+        m.print_infeasibilities()
         raise RuntimeError("Solving status 'infeasible'")
 
     return n