diff --git a/rules/build_sector.smk b/rules/build_sector.smk
index 1e8c70ba..ab8ff4ed 100644
--- a/rules/build_sector.smk
+++ b/rules/build_sector.smk
@@ -90,6 +90,7 @@ if config["sector"]["gas_network"] or config["sector"]["H2_retrofit"]:
                 keep_local=True,
             ),
             entry="data/gas_network/scigrid-gas/data/IGGIELGN_BorderPoints.geojson",
+            storage="data/gas_network/scigrid-gas/data/IGGIELGN_Storages.geojson",
             regions_onshore=RESOURCES
             + "regions_onshore_elec_s{simpl}_{clusters}.geojson",
             regions_offshore=RESOURCES
diff --git a/scripts/build_gas_input_locations.py b/scripts/build_gas_input_locations.py
index 07707658..ad449202 100644
--- a/scripts/build_gas_input_locations.py
+++ b/scripts/build_gas_input_locations.py
@@ -75,7 +75,7 @@ def build_gem_prod_data(fn):
     return gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
 
 
-def build_gas_input_locations(gem_fn, entry_fn, countries):
+def build_gas_input_locations(gem_fn, entry_fn, sto_fn, countries):
     # LNG terminals
     lng = build_gem_lng_data(gem_fn)
 
@@ -88,23 +88,30 @@ def build_gas_input_locations(gem_fn, entry_fn, countries):
         | (entry.from_country == "NO")  # malformed datapoint  # entries from NO to GB
     ]
 
+    sto =  read_scigrid_gas(sto_fn)
+    remove_country = ["RU", "UA", "TR", "BY"]
+    sto = sto.query("country_code != @remove_country")
+
     # production sites inside the model scope
     prod = build_gem_prod_data(gem_fn)
 
     mcm_per_day_to_mw = 437.5  # MCM/day to MWh/h
     mcm_per_year_to_mw = 1.199  #  MCM/year to MWh/h
     mtpa_to_mw = 1649.224  # mtpa to MWh/h
-    lng["p_nom"] = lng["CapacityInMtpa"] * mtpa_to_mw
-    entry["p_nom"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw
-    prod["p_nom"] = prod["mcm_per_year"] * mcm_per_year_to_mw
+    mcm_to_gwh = 11.36  # MCM to GWh
+    lng["capacity"] = lng["CapacityInMtpa"] * mtpa_to_mw
+    entry["capacity"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw
+    prod["capacity"] = prod["mcm_per_year"] * mcm_per_year_to_mw
+    sto["capacity"] = sto["max_cushionGas_M_m3"] * mcm_to_gwh
 
     lng["type"] = "lng"
     entry["type"] = "pipeline"
     prod["type"] = "production"
+    sto["type"] = "storage"
 
-    sel = ["geometry", "p_nom", "type"]
+    sel = ["geometry", "capacity", "type"]
 
-    return pd.concat([prod[sel], entry[sel], lng[sel]], ignore_index=True)
+    return pd.concat([prod[sel], entry[sel], lng[sel], sto[sel]], ignore_index=True)
 
 
 if __name__ == "__main__":
@@ -137,6 +144,7 @@ if __name__ == "__main__":
     gas_input_locations = build_gas_input_locations(
         snakemake.input.gem,
         snakemake.input.entry,
+        snakemake.input.storage,
         countries,
     )
 
@@ -147,8 +155,8 @@ if __name__ == "__main__":
     gas_input_nodes.to_file(snakemake.output.gas_input_nodes, driver="GeoJSON")
 
     gas_input_nodes_s = (
-        gas_input_nodes.groupby(["bus", "type"])["p_nom"].sum().unstack()
+        gas_input_nodes.groupby(["bus", "type"])["capacity"].sum().unstack()
     )
-    gas_input_nodes_s.columns.name = "p_nom"
+    gas_input_nodes_s.columns.name = "capacity"
 
     gas_input_nodes_s.to_csv(snakemake.output.gas_input_nodes_simplified)
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index e9d97ade..9387d4b1 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -454,10 +454,11 @@ def add_carrier_buses(n, carrier, nodes=None):
     n.add("Carrier", carrier)
 
     unit = "MWh_LHV" if carrier == "gas" else "MWh_th"
+    # preliminary value for non-gas carriers to avoid zeros
+    capital_cost = costs.at["gas storage", "fixed"] if carrier == "gas" else 0.02
 
     n.madd("Bus", nodes, location=location, carrier=carrier, unit=unit)
 
-    # capital cost could be corrected to e.g. 0.2 EUR/kWh * annuity and O&M
     n.madd(
         "Store",
         nodes + " Store",
@@ -465,8 +466,7 @@ def add_carrier_buses(n, carrier, nodes=None):
         e_nom_extendable=True,
         e_cyclic=True,
         carrier=carrier,
-        capital_cost=0.2
-        * costs.at[carrier, "discount rate"],  # preliminary value to avoid zeros
+        capital_cost=capital_cost,
     )
 
     n.madd(
@@ -1162,7 +1162,7 @@ def add_storage_and_grids(n, costs):
 
     if options["gas_network"]:
         logger.info(
-            "Add natural gas infrastructure, incl. LNG terminals, production and entry-points."
+            "Add natural gas infrastructure, incl. LNG terminals, production, storage and entry-points."
         )
 
         if options["H2_retrofit"]:
@@ -1207,10 +1207,17 @@ def add_storage_and_grids(n, costs):
         remove_i = n.generators[gas_i & internal_i].index
         n.generators.drop(remove_i, inplace=True)
 
-        p_nom = gas_input_nodes.sum(axis=1).rename(lambda x: x + " gas")
+        input_types = ["lng", "pipeline", "production"]
+        p_nom = gas_input_nodes[input_types].sum(axis=1).rename(lambda x: x + " gas")
         n.generators.loc[gas_i, "p_nom_extendable"] = False
         n.generators.loc[gas_i, "p_nom"] = p_nom
 
+        # add existing gas storage capacity
+        gas_i = n.stores.carrier == "gas"
+        e_nom = gas_input_nodes["storage"].rename(lambda x: x + " gas Store").reindex(n.stores.index).fillna(0.) * 1e3 # MWh_LHV
+        e_nom.clip(upper=e_nom.quantile(0.98), inplace=True) # limit extremely large storage
+        n.stores.loc[gas_i, "e_nom_min"] = e_nom
+
         # add candidates for new gas pipelines to achieve full connectivity
 
         G = nx.Graph()