diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8dfc2f33f..fd39e428e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -25,6 +25,25 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+## Develop LDES
+### Added
+- Added input option **can_export_to_grid** (defaults to _false_) to `ElectricStorage` and decision variable **dvStorageToGrid**
+- Added input option **optimize_soc_init_fraction** (defaults to _false_) to `ElectricStorage`, which makes the optimization choose the inital SOC (equal to final SOC) instead of using **soc_init_fraction**. The initial SOC is also constrained to equal the final SOC, which eliminates the "free energy" issue. We currently do not fix SOC when **soc_init_fraction** is used because this has caused infeasibility. 
+- Added output **initial_capital_cost** to all techs
+- Added input **soc_based_per_ts_thermal_decay_fraction** to both `ColdThermalStorage` and `HotThermalStorage` for modeling SOC-based thermal decay rate 
+- Renamed **thermal_decay_rate_fraction** to **capacity_based_per_ts_thermal_decay_fraction** in both `ColdThermalStorage` and `HotThermalStorage` to clarify that the variable represents the capacity-based thermal decay rate
+- Added ability to model multiple `ElectricStorage` at once. Inputs are provided as a _list_ of _dict_ under the `ElectricStorage` key and outputs are similarly formatted. This mirrors the existing multiple `PV` implementation. This functionality is not yet implemented for MPC scenarios.
+### Fixed
+- Added missing outputs **lifecycle_export_benefit_before_tax** and **year_one_export_benefit_after_tax** to `ElectricTariff`
+- Add missing output **year_one_om_cost_before_tax** to `PV`
+
+## Develop battery_om
+### Added
+- Added inputs **om_cost_per_kw** and **om_cost_per_kwh** to `ElectricStorage` for modeling capacity-based O&M 
+- Added outputs **lifecycle_om_cost_after_tax** and **year_one_om_cost_before_tax** to `ElectricStorage` 
+- Added input **soc_based_per_ts_self_discharge_fraction** to `ElectricStorage` for modeling SOC-based battery self-discharge
+- Added input **capacity_based_per_ts_self_discharge_fraction** to `ElectricStorage` for modeling capacity-based battery self-discharge
+
 ## Develop - 2024-10-11
 ### Added
 - Add new **ElectricStorage** parameters **max_duration_hours** and **min_duration_hours** to bound the energy duration of battery storage
diff --git a/src/constraints/battery_degradation.jl b/src/constraints/battery_degradation.jl
index 4b31e5297..7f81a8326 100644
--- a/src/constraints/battery_degradation.jl
+++ b/src/constraints/battery_degradation.jl
@@ -11,7 +11,7 @@ function add_degradation_variables(m, p)
 end
 
 
-function constrain_degradation_variables(m, p; b="ElectricStorage")
+function constrain_degradation_variables(m, p, b)
     days = 1:365*p.s.financial.analysis_years
     ts_per_day = 24 / p.hours_per_time_step
     ts_per_year = ts_per_day * 365
@@ -44,14 +44,12 @@ end
 
 
 """
-    add_degradation(m, p; b="ElectricStorage")
+    add_degradation(m, p, b)
 
 NOTE the average SOC and EFC variables are in absolute units. For example, the SOH variable starts 
     at the battery capacity in kWh.
 """
-function add_degradation(m, p; b="ElectricStorage")
-    
-    # Indices
+function add_degradation(m, p, b)
     days = 1:365*p.s.financial.analysis_years
     months = 1:p.s.financial.analysis_years*12
 
@@ -75,7 +73,7 @@ function add_degradation(m, p; b="ElectricStorage")
     end
 
     add_degradation_variables(m, p)
-    constrain_degradation_variables(m, p, b=b)
+    constrain_degradation_variables(m, p, b)
 
     @constraint(m, [d in 2:days[end]],
         m[:SOH][d] == m[:SOH][d-1] - p.hours_per_time_step * (
diff --git a/src/constraints/electric_utility_constraints.jl b/src/constraints/electric_utility_constraints.jl
index c1842f578..d60a0f22d 100644
--- a/src/constraints/electric_utility_constraints.jl
+++ b/src/constraints/electric_utility_constraints.jl
@@ -129,8 +129,14 @@ function add_export_constraints(m, p; _n="")
         if typeof(binNEM) <: Real  # no need for wholesale binary
             binWHL = 1
             WHL_benefit = @expression(m, p.pwf_e * p.hours_per_time_step *
-                sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
-                        for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                #sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
+                #        for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                sum(p.s.electric_tariff.export_rates[:WHL][ts] *
+                    (
+                        sum(m[Symbol("dvStorageToGrid"*_n)][b, ts] for b in p.s.storage.types.elec) + 
+                        sum(m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] for t in p.techs_by_exportbin[:WHL])
+                    ) for ts in p.time_steps
+                )
             )
         else
             binWHL = @variable(m, binary = true)
@@ -273,13 +279,13 @@ function add_simultaneous_export_import_constraint(m, p; _n="")
             }
         )
     else
-        bigM_hourly_load = maximum(p.s.electric_load.loads_kw)+maximum(p.s.space_heating_load.loads_kw)+maximum(p.s.process_heat_load.loads_kw)+maximum(p.s.dhw_load.loads_kw)+maximum(p.s.cooling_load.loads_kw_thermal)
+        bigM_hourly_load_plus_battery_power = maximum(p.s.electric_load.loads_kw)+maximum(p.s.space_heating_load.loads_kw)+maximum(p.s.process_heat_load.loads_kw)+maximum(p.s.dhw_load.loads_kw)+maximum(p.s.cooling_load.loads_kw_thermal)+sum(Real[p.s.storage.attr[b].max_kw for b in p.s.storage.types.elec])
         @constraint(m, NoGridPurchasesBinary[ts in p.time_steps],
             sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) +
-            sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) <= bigM_hourly_load*(1-m[Symbol("binNoGridPurchases"*_n)][ts])
+            sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) <= bigM_hourly_load_plus_battery_power*(1-m[Symbol("binNoGridPurchases"*_n)][ts])
         )
         @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
-            sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) <= bigM_hourly_load * m[Symbol("binNoGridPurchases"*_n)][ts]
+            sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) <= bigM_hourly_load_plus_battery_power * m[Symbol("binNoGridPurchases"*_n)][ts]
         )
     end
 end
diff --git a/src/constraints/outage_constraints.jl b/src/constraints/outage_constraints.jl
index 2fdc6e92a..d9418df53 100644
--- a/src/constraints/outage_constraints.jl
+++ b/src/constraints/outage_constraints.jl
@@ -273,56 +273,56 @@ function add_binMGCHPIsOnInTS_constraints(m, p; _n="")
     # TODO? make binMGCHPIsOnInTS indexed on p.techs.chp    
 end
 
-function add_MG_storage_dispatch_constraints(m,p)
+function add_MG_elec_storage_dispatch_constraints(m,p,b)
     # initial SOC at start of each outage equals the grid-optimal SOC
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps],
-        m[:dvMGStoredEnergy][s, tz, 0] <= m[:dvStoredEnergy]["ElectricStorage", tz]
+        m[:dvMGStoredEnergy][s, tz, 0] <= m[:dvStoredEnergy][b, tz]
     )
     
     # state of charge
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
         m[:dvMGStoredEnergy][s, tz, ts] == m[:dvMGStoredEnergy][s, tz, ts-1] + p.hours_per_time_step * (
-            p.s.storage.attr["ElectricStorage"].charge_efficiency * sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec)
-            - m[:dvMGDischargeFromStorage][s, tz, ts] / p.s.storage.attr["ElectricStorage"].discharge_efficiency
+            p.s.storage.attr[b].charge_efficiency * sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec)
+            - m[:dvMGDischargeFromStorage][s, tz, ts] / p.s.storage.attr[b].discharge_efficiency
         )
     )
 
 	# Prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
     @constraint(m, [ts in p.time_steps_without_grid],
-        m[:dvStorageEnergy]["ElectricStorage"] >= m[:dvMGStoredEnergy][s, tz, ts-1] + p.hours_per_time_step * (  
-            p.s.storage.attr["ElectricStorage"].charge_efficiency * sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec) 
+        m[:dvStorageEnergy][b] >= m[:dvMGStoredEnergy][s, tz, ts-1] + p.hours_per_time_step * (  
+            p.s.storage.attr[b].charge_efficiency * sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec) 
         )
     )
 
     # Min SOC
-    if p.s.storage.attr["ElectricStorage"].soc_min_applies_during_outages
+    if p.s.storage.attr[b].soc_min_applies_during_outages
         # Minimum state of charge
         @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-            m[:dvMGStoredEnergy][s, tz, ts] >=  p.s.storage.attr["ElectricStorage"].soc_min_fraction * m[:dvStorageEnergy]["ElectricStorage"]
+            m[:dvMGStoredEnergy][s, tz, ts] >=  p.s.storage.attr[b].soc_min_fraction * m[:dvStorageEnergy][b]
         )
     end
     
     # Dispatch to MG electrical storage is no greater than inverter capacity
     # and can't charge the battery unless binMGStorageUsed = 1
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-        m[:dvStoragePower]["ElectricStorage"] >= sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec)
+        m[:dvStoragePower][b] >= sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec)
     )
     
     # Dispatch from MG storage is no greater than inverter capacity
     # and can't discharge from storage unless binMGStorageUsed = 1
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-        m[:dvStoragePower]["ElectricStorage"] >= m[:dvMGDischargeFromStorage][s, tz, ts]
+        m[:dvStoragePower][b] >= m[:dvMGDischargeFromStorage][s, tz, ts]
     )
     
     # Dispatch to and from electrical storage is no greater than power capacity
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-        m[:dvStoragePower]["ElectricStorage"] >= m[:dvMGDischargeFromStorage][s, tz, ts]
+        m[:dvStoragePower][b] >= m[:dvMGDischargeFromStorage][s, tz, ts]
             + sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec)
     )
     
     # State of charge upper bound is storage system size
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-        m[:dvStorageEnergy]["ElectricStorage"] >= m[:dvMGStoredEnergy][s, tz, ts]
+        m[:dvStorageEnergy][b] >= m[:dvMGStoredEnergy][s, tz, ts]
     )
     
     if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
@@ -335,11 +335,11 @@ function add_MG_storage_dispatch_constraints(m,p)
         )
     else
         @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-            sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec) <= p.s.storage.attr["ElectricStorage"].max_kw * m[:binMGStorageUsed]
+            sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec) <= p.s.storage.attr[b].max_kw * m[:binMGStorageUsed]
         )
         
         @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-            m[:dvMGDischargeFromStorage][s, tz, ts] <= p.s.storage.attr["ElectricStorage"].max_kw * m[:binMGStorageUsed]
+            m[:dvMGDischargeFromStorage][s, tz, ts] <= p.s.storage.attr[b].max_kw * m[:binMGStorageUsed]
         )
     end
 end
diff --git a/src/constraints/storage_constraints.jl b/src/constraints/storage_constraints.jl
index c46d4e775..9349217ef 100644
--- a/src/constraints/storage_constraints.jl
+++ b/src/constraints/storage_constraints.jl
@@ -27,21 +27,29 @@ function add_storage_size_constraints(m, p, b; _n="")
         @constraint(m,
             m[Symbol("dvStorageEnergy"*_n)][b] <= m[Symbol("dvStoragePower"*_n)][b] * p.s.storage.attr[b].max_duration_hours
         )
-
         @constraint(m,
             m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoragePower"*_n)][b] * p.s.storage.attr[b].min_duration_hours
         )
-    end    
-
-
+    end
 end
 
 
 function add_general_storage_dispatch_constraints(m, p, b; _n="")
-    # Constraint (4a): initial state of charge
-	@constraint(m,
-        m[Symbol("dvStoredEnergy"*_n)][b, 0] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
-    )
+    # Constraint (4a): initial and final state of charge
+    if hasproperty(p.s.storage.attr[b], :optimize_soc_init_fraction) && p.s.storage.attr[b].optimize_soc_init_fraction
+        # @constraint(m, m[:dvStoredEnergy][b,maximum(p.time_steps)] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b] )
+        print("\nOptimizing "*b*" inital SOC and constraining initial SOC = final SOC\n")
+        @constraint(m,
+            m[Symbol("dvStoredEnergy"*_n)][b, 0] == m[:dvStoredEnergy][b, maximum(p.time_steps)]
+        )
+    else
+        @constraint(m,
+            m[Symbol("dvStoredEnergy"*_n)][b, 0] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+        )
+        # @constraint(m,
+        #     m[Symbol("dvStoredEnergy"*_n)][b, maximum(p.time_steps)] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+        # )
+    end
 
     #Constraint (4n): State of charge upper bound is storage system size
     @constraint(m, [ts in p.time_steps],
@@ -65,27 +73,36 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
 				
 	# Constraint (4g)-1: state-of-charge for electrical storage - with grid
 	@constraint(m, [ts in p.time_steps_with_grid],
-        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
+        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + 
+        p.hours_per_time_step * (  
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) 
             + p.s.storage.attr[b].grid_charge_efficiency * m[Symbol("dvGridToStorage"*_n)][b, ts] 
-            - m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.s.storage.attr[b].discharge_efficiency
+            - ((m[Symbol("dvDischargeFromStorage"*_n)][b,ts]+m[Symbol("dvStorageToGrid"*_n)][b, ts])  / p.s.storage.attr[b].discharge_efficiency)
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_self_discharge_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_self_discharge_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
 	)
 	# Constraint (4g)-2: state-of-charge for electrical storage - no grid
 	@constraint(m, [ts in p.time_steps_without_grid],
-        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
+        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + 
+        p.hours_per_time_step * (  
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.elec) 
             - m[Symbol("dvDischargeFromStorage"*_n)][b, ts] / p.s.storage.attr[b].discharge_efficiency
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_self_discharge_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_self_discharge_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
     )	
 
 	# Constraint (4h): prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
     # (4h)-1: with grid
 	@constraint(m, [ts in p.time_steps_with_grid],
-        m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
+        m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + 
+        p.hours_per_time_step * (  
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) 
             + p.s.storage.attr[b].grid_charge_efficiency * m[Symbol("dvGridToStorage"*_n)][b, ts] 
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_self_discharge_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_self_discharge_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
 	)
 	# (4h)-2: no grid
 	@constraint(m, [ts in p.time_steps_without_grid],
@@ -104,21 +121,33 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
 	@constraint(m, [ts in p.time_steps_with_grid],
         m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts] + 
             sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) + m[Symbol("dvGridToStorage"*_n)][b, ts]
+            + m[Symbol("dvStorageToGrid"*_n)][b, ts]
     )
 
+    #Dispatch from electrical storage is no greater than power capacity 
+    @constraint(m, [ts in p.time_steps_without_grid],
+        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + m[Symbol("dvStorageToGrid"*_n)][b, ts])
+
 	#Constraint (4l)-alt: Dispatch from electrical storage is no greater than power capacity (no grid connection)
 	@constraint(m, [ts in p.time_steps_without_grid],
         m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + 
             sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec)
     )
-					
-    # Remove grid-to-storage as an option if option to grid charge is turned off
+			
+    #Constraint (4m)-1: Remove grid-to-storage as an option if option to grid charge is turned off
     if !(p.s.storage.attr[b].can_grid_charge)
         for ts in p.time_steps_with_grid
             fix(m[Symbol("dvGridToStorage"*_n)][b, ts], 0.0, force=true)
         end
 	end
 
+    #Constraint (4m)-2: Force storage export to grid to zero if option to grid export is turned off
+    if !p.s.storage.attr[b].can_export_to_grid
+        for ts in p.time_steps
+            fix(m[Symbol("dvStorageToGrid"*_n)][b, ts], 0.0, force=true)
+        end
+    end
+
     if p.s.storage.attr[b].minimum_avg_soc_fraction > 0
         avg_soc = sum(m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps) /
                    (8760. / p.hours_per_time_step)
@@ -134,17 +163,19 @@ function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
 	@constraint(m, [ts in p.time_steps],
         m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (
             p.s.storage.attr[b].charge_efficiency * sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads) -
-            sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] for q in p.heating_loads) / p.s.storage.attr[b].discharge_efficiency -
-            p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+            sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] for q in p.heating_loads) / p.s.storage.attr[b].discharge_efficiency
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_thermal_decay_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_thermal_decay_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
     )
     
     # Prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
 	@constraint(m, [ts in p.time_steps],
         m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (  
             p.s.storage.attr[b].charge_efficiency * sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads)
-            - p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_thermal_decay_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_thermal_decay_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
     )
 
     #Constraint (4n)-1: Dispatch to and from thermal storage is no greater than power capacity
@@ -185,17 +216,19 @@ function add_cold_thermal_storage_dispatch_constraints(m, p, b; _n="")
 	@constraint(m, ColdTESInventoryCon[ts in p.time_steps],
         m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling) -
-            m[Symbol("dvDischargeFromStorage"*_n)][b,ts]/p.s.storage.attr[b].discharge_efficiency -
-            p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+            m[Symbol("dvDischargeFromStorage"*_n)][b,ts]/p.s.storage.attr[b].discharge_efficiency
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_thermal_decay_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_thermal_decay_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
     )
 
     # Prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
 	@constraint(m, [ts in p.time_steps],
         m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (  
             p.s.storage.attr[b].charge_efficiency * sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling)
-            - p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
         )
+        - (p.s.storage.attr[b].soc_based_per_ts_thermal_decay_fraction * m[Symbol("dvStoredEnergy"*_n)][b, ts])
+        - (p.s.storage.attr[b].capacity_based_per_ts_thermal_decay_fraction * m[Symbol("dvStorageEnergy"*_n)][b])
     )
 
     #Constraint (4n)-2: Dispatch to and from thermal storage is no greater than power capacity
diff --git a/src/core/energy_storage/electric_storage.jl b/src/core/energy_storage/electric_storage.jl
index 771a84364..1579f3f21 100644
--- a/src/core/energy_storage/electric_storage.jl
+++ b/src/core/energy_storage/electric_storage.jl
@@ -167,10 +167,13 @@ end
     soc_min_applies_during_outages::Bool = false
     soc_init_fraction::Float64 = off_grid_flag ? 1.0 : 0.5
     can_grid_charge::Bool = off_grid_flag ? false : true
+    can_export_to_grid::Bool = false
     installed_cost_per_kw::Real = 910.0
     installed_cost_per_kwh::Real = 455.0
     replace_cost_per_kw::Real = 715.0
     replace_cost_per_kwh::Real = 318.0
+    om_cost_per_kw::Real = 0.0 # Capacity-based O&M costs in \$/kW-rated/year. When including battery replacement costs, the user should ensure that O&M costs do not account for augmentation/replacement.
+    om_cost_per_kwh::Real = 0.0 # Capacity-based O&M costs in \$/kWh-rated/year. When including battery replacement costs, the user should ensure that O&M costs do not account for augmentation/replacement.
     inverter_replacement_year::Int = 10
     battery_replacement_year::Int = 10
     macrs_option_years::Int = 7
@@ -185,6 +188,10 @@ end
     model_degradation::Bool = false
     degradation::Dict = Dict()
     minimum_avg_soc_fraction::Float64 = 0.0
+    capacity_based_per_ts_self_discharge_fraction::Float64 = 0.0 # Battery self-discharge as a fraction per timestep loss based on the rated kWh capacity of the sized storage system. 
+    soc_based_per_ts_self_discharge_fraction::Float64 = 0.0 # Battery self-discharge as a fraction per timestep loss based on kWh stored in each timestep
+    optimize_soc_init_fraction::Bool = false
+    name::String = "ElectricStorage"
     min_duration_hours::Real = 0.0 # Minimum amount of time storage can discharge at its rated power capacity
     max_duration_hours::Real = 100000.0 # Maximum amount of time storage can discharge at its rated power capacity (ratio of ElectricStorage size_kwh to size_kw)
 ```
@@ -202,10 +209,13 @@ Base.@kwdef struct ElectricStorageDefaults
     soc_min_applies_during_outages::Bool = false
     soc_init_fraction::Float64 = off_grid_flag ? 1.0 : 0.5
     can_grid_charge::Bool = off_grid_flag ? false : true
+    can_export_to_grid::Bool = false
     installed_cost_per_kw::Real = 910.0
     installed_cost_per_kwh::Real = 455.0
     replace_cost_per_kw::Real = 715.0
     replace_cost_per_kwh::Real = 318.0
+    om_cost_per_kw::Real = 0.0
+    om_cost_per_kwh::Real = 0.0
     inverter_replacement_year::Int = 10
     battery_replacement_year::Int = 10
     macrs_option_years::Int = 7
@@ -220,6 +230,10 @@ Base.@kwdef struct ElectricStorageDefaults
     model_degradation::Bool = false
     degradation::Dict = Dict()
     minimum_avg_soc_fraction::Float64 = 0.0
+    capacity_based_per_ts_self_discharge_fraction::Float64 = 0.0
+    soc_based_per_ts_self_discharge_fraction::Float64 = 0.0
+    optimize_soc_init_fraction::Bool = false
+    name::String = "ElectricStorage"
     min_duration_hours::Real = 0.0
     max_duration_hours::Real = 100000.0
 end
@@ -243,10 +257,13 @@ struct ElectricStorage <: AbstractElectricStorage
     soc_min_applies_during_outages::Bool
     soc_init_fraction::Float64
     can_grid_charge::Bool
+    can_export_to_grid::Bool
     installed_cost_per_kw::Real
     installed_cost_per_kwh::Real
     replace_cost_per_kw::Real
     replace_cost_per_kwh::Real
+    om_cost_per_kw::Real
+    om_cost_per_kwh::Real
     inverter_replacement_year::Int
     battery_replacement_year::Int
     macrs_option_years::Int
@@ -263,6 +280,10 @@ struct ElectricStorage <: AbstractElectricStorage
     model_degradation::Bool
     degradation::Degradation
     minimum_avg_soc_fraction::Float64
+    capacity_based_per_ts_self_discharge_fraction::Float64
+    soc_based_per_ts_self_discharge_fraction::Float64
+    optimize_soc_init_fraction::Bool
+    name::String
     min_duration_hours::Real
     max_duration_hours::Real
 
@@ -336,10 +357,13 @@ struct ElectricStorage <: AbstractElectricStorage
             s.soc_min_applies_during_outages,
             s.soc_init_fraction,
             s.can_grid_charge,
+            s.can_export_to_grid,
             s.installed_cost_per_kw,
             s.installed_cost_per_kwh,
             replace_cost_per_kw,
             replace_cost_per_kwh,
+            s.om_cost_per_kw,
+            s.om_cost_per_kwh,
             s.inverter_replacement_year,
             s.battery_replacement_year,
             s.macrs_option_years,
@@ -356,6 +380,10 @@ struct ElectricStorage <: AbstractElectricStorage
             s.model_degradation,
             degr,
             s.minimum_avg_soc_fraction,
+            s.capacity_based_per_ts_self_discharge_fraction,
+            s.soc_based_per_ts_self_discharge_fraction,
+            s.optimize_soc_init_fraction,
+            s.name,
             s.min_duration_hours,
             s.max_duration_hours
         )
diff --git a/src/core/energy_storage/thermal_storage.jl b/src/core/energy_storage/thermal_storage.jl
index 41ffb6aff..a0703578f 100644
--- a/src/core/energy_storage/thermal_storage.jl
+++ b/src/core/energy_storage/thermal_storage.jl
@@ -15,7 +15,8 @@ Cold thermal energy storage sytem; specifically, a chilled water system used to
     soc_min_fraction::Float64 = 0.1 # Minimum allowable TES thermal state of charge
     soc_init_fraction::Float64 = 0.5 # TES thermal state of charge at first hour of optimization
     installed_cost_per_gal::Float64 = 1.50 # Thermal energy-based cost of TES (e.g. volume of the tank)
-    thermal_decay_rate_fraction::Float64 = 0.0004 # Thermal loss (gain) rate as a fraction of energy storage capacity, per hour (frac*energy_capacity/hr = kw_thermal)
+    soc_based_per_ts_thermal_decay_fraction = 0.0 # Thermal loss (gain) rate as a fraction of energy storage capacity per timestep (frac*energy_capacity = kw_thermal); the provided default is for an hourly thermal loss/gain
+    capacity_based_per_ts_thermal_decay_fraction = 0.0004 # Thermal loss (gain) rate as a fraction of the energy storaged in each timestep (frac*energy_stored = kw_thermal)
     om_cost_per_gal::Float64 = 0.0 # Yearly fixed O&M cost dependent on storage energy size
     macrs_option_years::Int = 7
     macrs_bonus_fraction::Float64 = 0.6
@@ -33,7 +34,8 @@ Base.@kwdef struct ColdThermalStorageDefaults <: AbstractThermalStorageDefaults
     soc_min_fraction::Float64 = 0.1
     soc_init_fraction::Float64 = 0.5
     installed_cost_per_gal::Float64 = 1.50
-    thermal_decay_rate_fraction::Float64 = 0.0004
+    soc_based_per_ts_thermal_decay_fraction::Float64 = 0.0
+    capacity_based_per_ts_thermal_decay_fraction::Float64 = 0.0004
     om_cost_per_gal::Float64 = 0.0
     macrs_option_years::Int = 7
     macrs_bonus_fraction::Float64 = 0.6
@@ -55,7 +57,8 @@ end
     soc_min_fraction::Float64 = 0.1
     soc_init_fraction::Float64 = 0.5
     installed_cost_per_gal::Float64 = 1.50
-    thermal_decay_rate_fraction::Float64 = 0.0004
+    soc_based_per_ts_thermal_decay_fraction = 0.0
+    capacity_based_per_ts_thermal_decay_fraction = 0.0004
     om_cost_per_gal::Float64 = 0.0
     macrs_option_years::Int = 7
     macrs_bonus_fraction::Float64 = 0.6
@@ -76,7 +79,8 @@ Base.@kwdef struct HotThermalStorageDefaults <: AbstractThermalStorageDefaults
     soc_min_fraction::Float64 = 0.1
     soc_init_fraction::Float64 = 0.5
     installed_cost_per_gal::Float64 = 1.50
-    thermal_decay_rate_fraction::Float64 = 0.0004
+    soc_based_per_ts_thermal_decay_fraction::Float64 = 0.0
+    capacity_based_per_ts_thermal_decay_fraction::Float64 = 0.0004
     om_cost_per_gal::Float64 = 0.0
     macrs_option_years::Int = 7
     macrs_bonus_fraction::Float64 = 0.6
@@ -105,7 +109,8 @@ struct ColdThermalStorage <: AbstractThermalStorage
     soc_min_fraction::Float64
     soc_init_fraction::Float64
     installed_cost_per_gal::Float64
-    thermal_decay_rate_fraction::Float64
+    soc_based_per_ts_thermal_decay_fraction::Float64
+    capacity_based_per_ts_thermal_decay_fraction::Float64
     om_cost_per_gal::Float64
     macrs_option_years::Int
     macrs_bonus_fraction::Float64
@@ -154,7 +159,8 @@ struct ColdThermalStorage <: AbstractThermalStorage
             s.soc_min_fraction,
             s.soc_init_fraction,
             s.installed_cost_per_gal,
-            s.thermal_decay_rate_fraction,
+            s.soc_based_per_ts_thermal_decay_fraction,
+            s.capacity_based_per_ts_thermal_decay_fraction,
             s.om_cost_per_gal,
             s.macrs_option_years,
             s.macrs_bonus_fraction,
@@ -188,7 +194,8 @@ struct HotThermalStorage <: AbstractThermalStorage
     soc_min_fraction::Float64
     soc_init_fraction::Float64
     installed_cost_per_gal::Float64
-    thermal_decay_rate_fraction::Float64
+    soc_based_per_ts_thermal_decay_fraction::Float64
+    capacity_based_per_ts_thermal_decay_fraction::Float64
     om_cost_per_gal::Float64
     macrs_option_years::Int
     macrs_bonus_fraction::Float64
@@ -240,7 +247,8 @@ struct HotThermalStorage <: AbstractThermalStorage
             s.soc_min_fraction,
             s.soc_init_fraction,
             s.installed_cost_per_gal,
-            s.thermal_decay_rate_fraction,
+            s.soc_based_per_ts_thermal_decay_fraction,
+            s.capacity_based_per_ts_thermal_decay_fraction,
             s.om_cost_per_gal,
             s.macrs_option_years,
             s.macrs_bonus_fraction,
diff --git a/src/core/reopt.jl b/src/core/reopt.jl
index 662d249cd..ca9e50317 100644
--- a/src/core/reopt.jl
+++ b/src/core/reopt.jl
@@ -151,6 +151,9 @@ function run_reopt(ms::AbstractArray{T, 1}, p::REoptInputs) where T <: JuMP.Abst
 			if !isempty(p.techs.pv)
 				organize_multiple_pv_results(p, results_dict)
 			end
+			if !isempty(p.s.storage.types.elec)
+				organize_multiple_elec_stor_results(p, results_dict)
+			end
 			return results_dict
 		else
 			throw(@error("REopt scenarios solved either with errors or non-optimal solutions."))
@@ -213,6 +216,7 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 				@constraint(m, [ts in p.time_steps], m[:dvGridToStorage][b, ts] == 0)
 				@constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
 						m[:dvProductionToStorage][b, t, ts] == 0)
+				@constraint(m, [ts in p.time_steps], m[Symbol("dvStorageToGrid"*_n)][b, ts] == 0)  # if there isn't a battery, then the battery can't export power to the grid
 			elseif b in p.s.storage.types.hot
 				@constraint(m, [q in q in setdiff(p.heating_loads, p.heating_loads_served_by_tes[b]), ts in p.time_steps], m[:dvHeatFromStorage][b,q,ts] == 0)
 				if "DomesticHotWater" in p.heating_loads_served_by_tes[b]
@@ -410,9 +414,11 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 		sum( p.s.storage.attr[b].net_present_cost_per_kwh * m[:dvStorageEnergy][b] for b in p.s.storage.types.all )
 	))
 	
-	@expression(m, TotalPerUnitSizeOMCosts, p.third_party_factor * p.pwf_om *
-		sum( p.om_cost_per_kw[t] * m[:dvSize][t] for t in p.techs.all )
-	)
+	@expression(m, TotalPerUnitSizeOMCosts, p.third_party_factor * p.pwf_om * (
+		sum(p.om_cost_per_kw[t] * m[:dvSize][t] for t in p.techs.all) + 
+		sum(p.s.storage.attr[b].om_cost_per_kw * m[:dvStoragePower][b] for b in p.s.storage.types.elec) +
+		sum(p.s.storage.attr[b].om_cost_per_kwh * m[:dvStorageEnergy][b] for b in p.s.storage.types.elec)
+	))
 
 	add_elec_utility_expressions(m, p)
 
@@ -420,11 +426,21 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
         add_dv_UnservedLoad_constraints(m,p)
 		add_outage_cost_constraints(m,p)
 		add_MG_production_constraints(m,p)
-		if !isempty(p.s.storage.types.elec)
-			add_MG_storage_dispatch_constraints(m,p)
+
+		if length(p.s.storage.types.elec) > 1
+			throw(@error("REopt does not support considering multiple ElectricStorage when modelling outages."))
+		end
+		if !isempty(p.s.storage.types.elec) #length=1
+			b = p.s.storage.types.elec[1]
+			if p.s.storage.attr[b].max_kw == 0 || p.s.storage.attr[b].max_kwh == 0
+				fix_MG_storage_variables(m,p)
+			else
+				add_MG_elec_storage_dispatch_constraints(m,p,b)
+			end
 		else
 			fix_MG_storage_variables(m,p)
 		end
+
 		add_cannot_have_MG_with_only_PVwind_constraints(m,p)
 		add_MG_size_constraints(m,p)
 		
@@ -470,7 +486,7 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 
 	for b in p.s.storage.types.elec
 		if p.s.storage.attr[b].model_degradation
-			add_degradation(m, p; b=b)
+			add_degradation(m, p, b)
 			if p.s.settings.add_soc_incentive # this warning should be tied to IF condition where SOC incentive is added
 				@warn "Settings.add_soc_incentive is set to true and it will incentivize BESS energy levels to be kept high. It could conflict with the battery degradation model and should be disabled."
 			end
@@ -586,6 +602,9 @@ function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true)
 		if organize_pvs && !isempty(p.techs.pv)  # do not want to organize_pvs when running BAU case in parallel b/c then proform code fails
 			organize_multiple_pv_results(p, results)
 		end
+		if !isempty(p.s.storage.types.elec)
+			organize_multiple_elec_stor_results(p, results)
+		end
 
 		# add error messages (if any) and warnings to results dict
 		results["Messages"] = logger_to_dict()
@@ -616,6 +635,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
 		dvProductionToStorage[p.s.storage.types.all, union(p.techs.ghp,p.techs.all), p.time_steps] >= 0  # Power from technology t used to charge storage system b [kW]
 		dvDischargeFromStorage[p.s.storage.types.all, p.time_steps] >= 0 # Power discharged from storage system b [kW]
 		dvGridToStorage[p.s.storage.types.elec, p.time_steps] >= 0 # Electrical power delivered to storage by the grid [kW]
+		dvStorageToGrid[p.s.storage.types.elec, p.time_steps] >= 0 # TODO, add: "p.StorageSalesTiers" as well? export of energy from storage to the grid
 		dvStoredEnergy[p.s.storage.types.all, 0:p.time_steps[end]] >= 0  # State of charge of storage system b
 		dvStoragePower[p.s.storage.types.all] >= 0   # Power capacity of storage system b [kW]
 		dvStorageEnergy[p.s.storage.types.all] >= 0   # Energy capacity of storage system b [kWh]
@@ -623,6 +643,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
 		dvPeakDemandMonth[p.months, 1:p.s.electric_tariff.n_monthly_demand_tiers] >= 0  # Peak electrical power demand during month m [kW]
 		MinChargeAdder >= 0
         binGHP[p.ghp_options], Bin  # Can be <= 1 if require_ghp_purchase=0, and is ==1 if require_ghp_purchase=1
+				
 	end
 
 	if !isempty(p.techs.gen)  # Problem becomes a MILP
@@ -641,7 +662,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
     end
 
 	if !(p.s.electric_utility.allow_simultaneous_export_import) & !isempty(p.s.electric_tariff.export_bins)
-		@warn "Adding binary variable to prevent simultaneous grid import/export. Some solvers are very slow with integer variables"
+		@warn "Adding binary variable to prevent simultaneous grid import/export. Some solvers are very slow with integer variables."
 		@variable(m, binNoGridPurchases[p.time_steps], Bin)
 	end
 
@@ -673,7 +694,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
     end
 
 	if !isempty(p.s.electric_utility.outage_durations) # add dvUnserved Load if there is at least one outage
-		@warn "Adding binary variable to model outages. Some solvers are very slow with integer variables"
+		@warn "Adding binary variable to model outages. Some solvers are very slow with integer variables."
 		max_outage_duration = maximum(p.s.electric_utility.outage_durations)
 		outage_time_steps = p.s.electric_utility.outage_time_steps
 		tZeros = p.s.electric_utility.outage_start_time_steps
diff --git a/src/core/reopt_multinode.jl b/src/core/reopt_multinode.jl
index a41a80a35..e49442693 100644
--- a/src/core/reopt_multinode.jl
+++ b/src/core/reopt_multinode.jl
@@ -2,7 +2,6 @@
 
 
 function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario
-	
 	dvs_idx_on_techs = String[
 		"dvSize",
 		"dvPurchaseSize",
@@ -16,7 +15,8 @@ function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}
 		"dvStorageEnergy",
 	]
 	dvs_idx_on_storagetypes_time_steps = String[
-		"dvDischargeFromStorage"
+		"dvDischargeFromStorage",
+		"dvStorageToGrid"
 	]
 	for p in ps
 		_n = string("_", p.s.site.node)
diff --git a/src/core/scenario.jl b/src/core/scenario.jl
index 621781f53..5fd8df026 100644
--- a/src/core/scenario.jl
+++ b/src/core/scenario.jl
@@ -61,7 +61,7 @@ A Scenario struct can contain the following keys:
 - [ASHPSpaceHeater](@ref) (optional)
 - [ASHPWaterHeater](@ref) (optional)
 
-All values of `d` are expected to be `Dicts` except for `PV` and `GHP`, which can be either a `Dict` or `Dict[]` (for multiple PV arrays or GHP options).
+All values of `d` are expected to be `Dicts` except for `PV`, `ElectricStorage`, and `GHP`, which can be either a `Dict` or `Dict[]` (for multiple PV arrays or GHP options).
 
 !!! note 
     Set `flex_hvac_from_json=true` if `FlexibleHVAC` values were loaded in from JSON (necessary to 
@@ -93,10 +93,11 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
     if haskey(d, "PV")
         if typeof(d["PV"]) <: AbstractArray
             for (i, pv) in enumerate(d["PV"])
-                if !(haskey(pv, "name"))
-                    pv["name"] = string("PV", i)
+                pv = dictkeys_tosymbols(pv)
+                if !(haskey(pv, :name))
+                    pv[:name] = string("PV", i)
                 end
-                push!(pvs, PV(;dictkeys_tosymbols(pv)..., off_grid_flag = settings.off_grid_flag, 
+                push!(pvs, PV(;pv..., off_grid_flag = settings.off_grid_flag, 
                             latitude=site.latitude))
             end
         elseif typeof(d["PV"]) <: AbstractDict
@@ -168,12 +169,27 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
         
     storage_structs = Dict{String, AbstractStorage}()
     if haskey(d,  "ElectricStorage")
-        storage_dict = dictkeys_tosymbols(d["ElectricStorage"])
-        storage_dict[:off_grid_flag] = settings.off_grid_flag
+        if typeof(d["ElectricStorage"]) <: AbstractArray
+            for (i, elecstor) in enumerate(d["ElectricStorage"])
+                storage_dict = dictkeys_tosymbols(elecstor)
+                storage_dict[:off_grid_flag] = settings.off_grid_flag
+                if !(haskey(storage_dict, :name))
+                    storage_dict[:name] = string("ElectricStorage", i)
+                end
+                storage_structs[storage_dict[:name]] = ElectricStorage(storage_dict, financial)
+            end
+        elseif typeof(d["ElectricStorage"]) <: AbstractDict
+            storage_dict = dictkeys_tosymbols(d["ElectricStorage"])
+            storage_dict[:off_grid_flag] = settings.off_grid_flag
+            storage_struct = ElectricStorage(storage_dict, financial)
+            storage_structs[storage_struct.name] = storage_struct
+        else
+            throw(@error("ElectricStorage input must be Dict or Dict[]."))
+        end
     else
         storage_dict = Dict(:max_kw => 0.0) 
+        storage_structs["ElectricStorage"] = ElectricStorage(storage_dict, financial)
     end
-    storage_structs["ElectricStorage"] = ElectricStorage(storage_dict, financial)
     # TODO stop building ElectricStorage when it is not modeled by user 
     #       (requires significant changes to constraints, variables)
     if haskey(d, "HotThermalStorage")
diff --git a/src/mpc/model.jl b/src/mpc/model.jl
index d40a80da3..43b992de2 100644
--- a/src/mpc/model.jl
+++ b/src/mpc/model.jl
@@ -34,10 +34,10 @@ Returns a Dict of results with keys matching those in the `MPCScenario`.
 function run_mpc(m::JuMP.AbstractModel, p::MPCInputs)
     build_mpc!(m, p)
 
-    if !p.s.settings.add_soc_incentive || !("ElectricStorage" in p.s.storage.types.elec)
+    if !p.s.settings.add_soc_incentive || isempty(p.s.storage.types.elec)
 		@objective(m, Min, m[:Costs])
 	else # Keep SOC high
-		@objective(m, Min, m[:Costs] - sum(m[:dvStoredEnergy]["ElectricStorage", ts] for ts in p.time_steps) /
+		@objective(m, Min, m[:Costs] - sum(m[:dvStoredEnergy][b, ts] for ts in p.time_steps, b in p.s.storage.types.elec) /
 									   (8760. / p.hours_per_time_step)
 		)
 	end
@@ -99,6 +99,7 @@ function build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)
 			@constraint(m, [ts in p.time_steps], m[:dvDischargeFromStorage][b, ts] == 0)
 			if b in p.s.storage.types.elec
 				@constraint(m, [ts in p.time_steps], m[:dvGridToStorage][b, ts] == 0)
+				@constraint(m, [ts in p.time_steps], m[:dvStorageToGrid][b, ts] == 0)
 			end
 		else
 			add_general_storage_dispatch_constraints(m, p, b)
@@ -170,16 +171,25 @@ function build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)
     add_previous_monthly_peak_constraint(m, p)
     add_previous_tou_peak_constraint(m, p)
 
-    # TODO: random outages in MPC?
 	if !isempty(p.s.electric_utility.outage_durations)
 		add_dv_UnservedLoad_constraints(m,p)
 		add_outage_cost_constraints(m,p)
 		add_MG_production_constraints(m,p)
-		if !isempty(p.s.storage.types.elec)	
-			add_MG_storage_dispatch_constraints(m,p)
+
+		if length(p.s.storage.types.elec) > 1
+			throw(@error("REopt does not support considering multiple ElectricStorage when modelling outages."))
+		end
+		if !isempty(p.s.storage.types.elec) #length=1
+			b = p.s.storage.types.elec[1]
+			if p.s.storage.attr[b].max_kw == 0 || p.s.storage.attr[b].max_kwh == 0
+				fix_MG_storage_variables(m,p)
+			else
+				add_MG_elec_storage_dispatch_constraints(m,p,b)
+			end
 		else
 			fix_MG_storage_variables(m,p)
 		end
+		
 		add_cannot_have_MG_with_only_PVwind_constraints(m,p)
 		add_MG_size_constraints(m,p)
 		
@@ -230,6 +240,7 @@ function add_variables!(m::JuMP.AbstractModel, p::MPCInputs)
 		dvCurtail[p.techs.all, p.time_steps] >= 0  # [kW]
 		dvProductionToStorage[p.s.storage.types.all, p.techs.all, p.time_steps] >= 0  # Power from technology t used to charge storage system b [kW]
 		dvDischargeFromStorage[p.s.storage.types.all, p.time_steps] >= 0 # Power discharged from storage system b [kW]
+		dvStorageToGrid[p.s.storage.types.elec, p.time_steps] >= 0 # TODO, add: "p.StorageSalesTiers" as well? export of energy from storage to the grid
 		dvGridToStorage[p.s.storage.types.elec, p.time_steps] >= 0 # Electrical power delivered to storage by the grid [kW]
 		dvStoredEnergy[p.s.storage.types.all, 0:p.time_steps[end]] >= 0  # State of charge of storage system b
 		dvStoragePower[p.s.storage.types.all] >= 0   # Power capacity of storage system b [kW]
@@ -248,8 +259,8 @@ function add_variables!(m::JuMP.AbstractModel, p::MPCInputs)
     m[:dvSize] = p.existing_sizes
 
 	for b in p.s.storage.types.all
-		fix(m[:dvStoragePower][b], p.s.storage.attr["ElectricStorage"].size_kw, force=true)
-		fix(m[:dvStorageEnergy][b], p.s.storage.attr["ElectricStorage"].size_kwh, force=true)
+		fix(m[:dvStoragePower][b], p.s.storage.attr[b].size_kw, force=true)
+		fix(m[:dvStorageEnergy][b], p.s.storage.attr[b].size_kwh, force=true)
 	end
 
 	# not modeling min charges since control does not affect them
diff --git a/src/mpc/model_multinode.jl b/src/mpc/model_multinode.jl
index 68a3ef34a..245dd3582 100644
--- a/src/mpc/model_multinode.jl
+++ b/src/mpc/model_multinode.jl
@@ -152,7 +152,8 @@ function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{MPCInputs})
 		"dvRatedProduction",
 	]
 	dvs_idx_on_storagetypes_time_steps = String[
-		"dvDischargeFromStorage"
+		"dvDischargeFromStorage",
+		"dvStorageToGrid"
 	]
 	for p in ps
 		_n = string("_", p.s.node)
diff --git a/src/mpc/structs.jl b/src/mpc/structs.jl
index 340bb07c0..47508e056 100644
--- a/src/mpc/structs.jl
+++ b/src/mpc/structs.jl
@@ -226,6 +226,7 @@ Base.@kwdef struct MPCElectricStorage < AbstractElectricStorage
     soc_min_fraction::Float64 = 0.2
     soc_init_fraction::Float64 = 0.5
     can_grid_charge::Bool = true
+    can_export_to_grid::Bool = false
     grid_charge_efficiency::Float64 = 0.96 * 0.975^2
 end
 ```
@@ -238,10 +239,13 @@ Base.@kwdef struct MPCElectricStorage <: AbstractElectricStorage
     soc_min_fraction::Float64 = 0.2
     soc_init_fraction::Float64 = 0.5
     can_grid_charge::Bool = true
+    can_export_to_grid::Bool = false
     grid_charge_efficiency::Float64 = 0.96 * 0.975^2
     max_kw::Float64 = size_kw
     max_kwh::Float64 = size_kwh
     minimum_avg_soc_fraction::Float64 = 0.0
+    capacity_based_per_ts_self_discharge_fraction::Float64 = 0.0
+    soc_based_per_ts_self_discharge_fraction::Float64 = 0.0
 end
 
 
diff --git a/src/outagesim/backup_reliability.jl b/src/outagesim/backup_reliability.jl
index d089d207c..1cd70e762 100644
--- a/src/outagesim/backup_reliability.jl
+++ b/src/outagesim/backup_reliability.jl
@@ -719,8 +719,14 @@ function backup_reliability_reopt_inputs(;d::Dict, p::REoptInputs, r::Dict = Dic
         !haskey(d, "Outages") ||
         Bool(get(d["Outages"], "electric_storage_microgrid_upgraded", false))
     )
-        r2[:battery_charge_efficiency] = p.s.storage.attr["ElectricStorage"].charge_efficiency
-        r2[:battery_discharge_efficiency] = p.s.storage.attr["ElectricStorage"].discharge_efficiency
+        if typeof(d["ElectricStorage"]) <: AbstractArray && length(d["ElectricStorage"]) == 1
+            d["ElectricStorage"] = d["ElectricStorage"][1]
+        elseif !(typeof(d["ElectricStorage"]) <: AbstractDict)
+            throw(@error("Calculating resilience metrics for a REopt solution with multiple ElectricStorage is not yet supported."))
+        end
+        elec_stor_name = get(d["ElectricStorage"],"name","ElectricStorage")
+        r2[:battery_charge_efficiency] = p.s.storage.attr[elec_stor_name].charge_efficiency
+        r2[:battery_discharge_efficiency] = p.s.storage.attr[elec_stor_name].discharge_efficiency
         r2[:battery_size_kw] = get(d["ElectricStorage"], "size_kw", 0)
 
         #ERP tool uses effective battery size so need to subtract minimum SOC
diff --git a/src/outagesim/outage_simulator.jl b/src/outagesim/outage_simulator.jl
index d365347dc..47a4d9826 100644
--- a/src/outagesim/outage_simulator.jl
+++ b/src/outagesim/outage_simulator.jl
@@ -223,9 +223,6 @@ Returns a dict
 ```
 """
 function simulate_outages(d::Dict, p::REoptInputs; microgrid_only::Bool=false)
-    batt_roundtrip_efficiency = (p.s.storage.attr["ElectricStorage"].charge_efficiency *
-                                p.s.storage.attr["ElectricStorage"].discharge_efficiency)
-
     # TODO handle generic PV names
     pv_kw_ac_hourly = zeros(length(p.time_steps))
     if "PV" in keys(d) && !(microgrid_only && !Bool(get(d["Outages"], "PV_upgraded", false)))
@@ -250,11 +247,20 @@ function simulate_outages(d::Dict, p::REoptInputs; microgrid_only::Bool=false)
     batt_kwh = 0
     batt_kw = 0
     init_soc = zeros(length(p.time_steps))
+    elec_stor_name = "ElectricStorage"
     if "ElectricStorage" in keys(d)
+        if typeof(d["ElectricStorage"]) <: AbstractArray && length(d["ElectricStorage"]) == 1
+            d["ElectricStorage"] = d["ElectricStorage"][1]
+        elseif !(typeof(d["ElectricStorage"]) <: AbstractDict)
+            throw(@error("simulate_outages for a REopt solution with multiple ElectricStorage is not yet supported"))
+        end
+        elec_stor_name = get(d["ElectricStorage"],"name","ElectricStorage")
         batt_kwh = get(d["ElectricStorage"], "size_kwh", 0)
         batt_kw = get(d["ElectricStorage"], "size_kw", 0)
         init_soc = get(d["ElectricStorage"], "soc_series_fraction", zeros(length(p.time_steps)))
     end
+    batt_roundtrip_efficiency = (p.s.storage.attr[elec_stor_name].charge_efficiency *
+                                p.s.storage.attr[elec_stor_name].discharge_efficiency)
     if microgrid_only && !Bool(get(d["Outages"], "electric_storage_microgrid_upgraded", false))
         batt_kwh = 0
         batt_kw = 0
diff --git a/src/results/absorption_chiller.jl b/src/results/absorption_chiller.jl
index 8bdd47069..af51ee435 100644
--- a/src/results/absorption_chiller.jl
+++ b/src/results/absorption_chiller.jl
@@ -25,6 +25,7 @@ function add_absorption_chiller_results(m::JuMP.AbstractModel, p::REoptInputs, d
 
 	r["size_kw"] = value(sum(m[:dvSize][t] for t in p.techs.absorption_chiller))
 	r["size_ton"] = r["size_kw"] / KWH_THERMAL_PER_TONHOUR
+	r["initial_capital_cost"] = round(value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.absorption_chiller)) * p.s.absorption_chiller.installed_cost_per_kw, digits=3)
 	@expression(m, ABSORPCHLtoTESKW[ts in p.time_steps],
 		sum(m[:dvProductionToStorage][b,t,ts] for b in p.s.storage.types.cold, t in p.techs.absorption_chiller))
 	r["thermal_to_storage_series_ton"] = round.(value.(ABSORPCHLtoTESKW) ./ KWH_THERMAL_PER_TONHOUR, digits=5)
diff --git a/src/results/boiler.jl b/src/results/boiler.jl
index cacdfb695..b4cf002b9 100644
--- a/src/results/boiler.jl
+++ b/src/results/boiler.jl
@@ -22,6 +22,7 @@
 function add_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r = Dict{String, Any}()
     r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["Boiler"]) / KWH_PER_MMBTU, digits=3)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)]["Boiler"]) * p.s.boiler.installed_cost_per_kw, digits=3)
 	r["fuel_consumption_series_mmbtu_per_hour"] = 
         round.(value.(m[:dvFuelUsage]["Boiler", ts] for ts in p.time_steps) / KWH_PER_MMBTU, digits=3)
     r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]), digits=3)
diff --git a/src/results/chp.jl b/src/results/chp.jl
index ccd5aa5f0..b154becbf 100644
--- a/src/results/chp.jl
+++ b/src/results/chp.jl
@@ -30,7 +30,8 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 	# Note: the node number is an empty string if evaluating a single `Site`.
     r = Dict{String, Any}()
 	r["size_kw"] = value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.chp))
-    r["size_supplemental_firing_kw"] = value(sum(m[Symbol("dvSupplementaryFiringSize"*_n)][t] for t in p.techs.chp))
+	#TODO: add initial_capital_cost to results (need to handle installed_cost_per_kw being a vector when cost curve is used)
+	r["size_supplemental_firing_kw"] = value(sum(m[Symbol("dvSupplementaryFiringSize"*_n)][t] for t in p.techs.chp))
 	@expression(m, CHPFuelUsedKWH, sum(m[Symbol("dvFuelUsage"*_n)][t, ts] for t in p.techs.chp, ts in p.time_steps))
 	r["annual_fuel_consumption_mmbtu"] = round(value(CHPFuelUsedKWH) / KWH_PER_MMBTU, digits=3)
 	@expression(m, Year1CHPElecProd,
@@ -59,7 +60,7 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r["electric_to_grid_series_kw"] = round.(value.(CHPtoGrid), digits=3)
 	if !isempty(p.s.storage.types.elec)
 		@expression(m, CHPtoBatt[ts in p.time_steps],
-			sum(m[Symbol("dvProductionToStorage"*_n)]["ElectricStorage",t,ts] for t in p.techs.chp))
+			sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.chp, b in p.s.storage.types.elec))
 	else
 		CHPtoBatt = zeros(length(p.time_steps))
 	end
diff --git a/src/results/electric_heater.jl b/src/results/electric_heater.jl
index 024804423..e537ce870 100644
--- a/src/results/electric_heater.jl
+++ b/src/results/electric_heater.jl
@@ -20,6 +20,7 @@
 function add_electric_heater_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r = Dict{String, Any}()
     r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["ElectricHeater"]) / KWH_PER_MMBTU, digits=3)
+    r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)]["ElectricHeater"]) * p.s.electric_heater.installed_cost_per_kw, digits=3)
     @expression(m, ElectricHeaterElectricConsumptionSeries[ts in p.time_steps],
         p.hours_per_time_step * sum(m[:dvHeatingProduction][t,q,ts] / p.heating_cop[t][ts] 
         for q in p.heating_loads, t in p.techs.electric_heater))
diff --git a/src/results/electric_storage.jl b/src/results/electric_storage.jl
index ef07eb48c..a9bc7164e 100644
--- a/src/results/electric_storage.jl
+++ b/src/results/electric_storage.jl
@@ -34,8 +34,16 @@ function add_electric_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::
         r["initial_capital_cost"] = r["size_kwh"] * p.s.storage.attr[b].installed_cost_per_kwh +
             r["size_kw"] * p.s.storage.attr[b].installed_cost_per_kw
 
+        StoragePerUnitOMCosts = p.third_party_factor * p.pwf_om * (p.s.storage.attr[b].om_cost_per_kw * m[Symbol("dvStoragePower"*_n)][b] +
+                                                                 p.s.storage.attr[b].om_cost_per_kwh * m[Symbol("dvStorageEnergy"*_n)][b])
+
+        r["lifecycle_om_cost_after_tax"] = round(value(StoragePerUnitOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
+        r["lifecycle_om_cost_before_tax"] = round(value(StoragePerUnitOMCosts), digits=0)
+        r["year_one_om_cost_before_tax"] = round(value(StoragePerUnitOMCosts) / (p.pwf_om * p.third_party_factor), digits=0)
+        r["year_one_om_cost_after_tax"] = round(value(StoragePerUnitOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction) / (p.pwf_om * p.third_party_factor), digits=0)
+            
         if p.s.storage.attr[b].model_degradation
-            r["state_of_health"] = value.(m[:SOH]).data / value.(m[:dvStorageEnergy])["ElectricStorage"];
+            r["state_of_health"] = value.(m[:SOH]).data / value.(m[:dvStorageEnergy])[b];
             r["maintenance_cost"] = value(m[:degr_cost])
             if p.s.storage.attr[b].degradation.maintenance_strategy == "replacement"
                 r["replacement_month"] = round(Int, value(
@@ -44,9 +52,13 @@ function add_electric_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::
             end
             r["residual_value"] = value(m[:residual_value])
          end
+         # report the exported electricity from the battery:
+         r["storage_to_grid_series_kw"] = round.(value.(m[Symbol("dvStorageToGrid"*_n)][b, ts] for ts in p.time_steps), digits = 3)
+
     else
         r["soc_series_fraction"] = []
         r["storage_to_load_series_kw"] = []
+        r["storage_to_grid_series_kw"] = []
     end
 
     d[b] = r
@@ -68,4 +80,24 @@ function add_electric_storage_results(m::JuMP.AbstractModel, p::MPCInputs, d::Di
 
     d[b] = r
     nothing
+end
+
+"""
+    organize_multiple_elec_stor_results(p::REoptInputs, d::Dict)
+
+The last step in results processing: if more than one ElectricStorage was modeled then move their results from the top
+level keys (that use each ElectricStorage.name) to an array of results with "ElectricStorage" as the top key in the results dict `d`.
+"""
+function organize_multiple_elec_stor_results(p::REoptInputs, d::Dict)
+    if length(p.s.storage.types.elec) == 1 && p.s.storage.types.elec[1] == "ElectricStorage"
+        return nothing
+    end
+    stors = Dict[]
+    for storname in p.s.storage.types.elec
+        d[storname]["name"] = storname  # add name to results dict to distinguish each ElectricStorage
+        push!(stors, d[storname])
+        delete!(d, storname)
+    end
+    d["ElectricStorage"] = stors
+    nothing
 end
\ No newline at end of file
diff --git a/src/results/electric_tariff.jl b/src/results/electric_tariff.jl
index 8bc2c05a9..331c07a8e 100644
--- a/src/results/electric_tariff.jl
+++ b/src/results/electric_tariff.jl
@@ -36,7 +36,9 @@ function add_electric_tariff_results(m::JuMP.AbstractModel, p::REoptInputs, d::D
     r["year_one_min_charge_adder_before_tax"] = round(value(m[Symbol("MinChargeAdder"*_n)]) / p.pwf_e, digits=2)
                                 
     r["lifecycle_export_benefit_after_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=2)
+    r["lifecycle_export_benefit_before_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]), digits=2)
     r["year_one_export_benefit_before_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) / p.pwf_e, digits=0)
+    r["year_one_export_benefit_after_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) * (1 - p.s.financial.offtaker_tax_rate_fraction) / p.pwf_e, digits=0)
 
     r["lifecycle_coincident_peak_cost_after_tax"] = round(value(m[Symbol("TotalCPCharges"*_n)]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=2)
     r["year_one_coincident_peak_cost_before_tax"] = round(value(m[Symbol("TotalCPCharges"*_n)]) / p.pwf_e, digits=2)
diff --git a/src/results/electric_utility.jl b/src/results/electric_utility.jl
index d5b21f467..d5e60805c 100644
--- a/src/results/electric_utility.jl
+++ b/src/results/electric_utility.jl
@@ -94,7 +94,7 @@ function add_electric_utility_results(m::JuMP.AbstractModel, p::MPCInputs, d::Di
                                                          tier in 1:p.s.electric_tariff.n_energy_tiers)
     r["energy_supplied_kwh"] = round(value(Year1UtilityEnergy), digits=2)
 
-    if p.s.storage.attr["ElectricStorage"].size_kwh > 0
+    if any([p.s.storage.attr[b].size_kwh for b in p.s.storage.types.elec] .> 0)
         GridToBatt = @expression(m, [ts in p.time_steps], 
             sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) 
 		)
diff --git a/src/results/generator.jl b/src/results/generator.jl
index ec9c2efb7..16a00bae1 100644
--- a/src/results/generator.jl
+++ b/src/results/generator.jl
@@ -32,6 +32,7 @@ function add_generator_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _
 			for t in p.techs.gen, ts in p.time_steps)
 	)
 	r["size_kw"] = round(value(sum(m[:dvSize][t] for t in p.techs.gen)), digits=2)
+	r["initial_capital_cost"] = round(value(sum(m[:dvSize][t] for t in p.techs.gen)) * p.s.generator.installed_cost_per_kw, digits=3)
 	r["lifecycle_fixed_om_cost_after_tax"] = round(value(GenPerUnitSizeOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
 	r["lifecycle_variable_om_cost_after_tax"] = round(value(m[:TotalPerUnitProdOMCosts]) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
 	r["lifecycle_fuel_cost_after_tax"] = round(value(m[:TotalGenFuelCosts]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=2)
@@ -89,7 +90,7 @@ function add_generator_results(m::JuMP.AbstractModel, p::MPCInputs, d::Dict; _n=
 	r["variable_om_cost"] = round(value(m[:TotalPerUnitProdOMCosts]), digits=0)
 	r["fuel_cost"] = round(value(m[:TotalGenFuelCosts]), digits=2)
 
-    if p.s.storage.attr["ElectricStorage"].size_kw > 0
+	if any([p.s.storage.attr[b].size_kw for b in p.s.storage.types.elec] .> 0)
         generatorToBatt = @expression(m, [ts in p.time_steps],
             sum(m[:dvProductionToStorage][b, t, ts] for b in p.s.storage.types.elec, t in p.techs.gen))
         r["to_battery_series_kw"] = round.(value.(generatorToBatt), digits=3).data
diff --git a/src/results/proforma.jl b/src/results/proforma.jl
index a765573e3..8959b4e58 100644
--- a/src/results/proforma.jl
+++ b/src/results/proforma.jl
@@ -69,30 +69,38 @@ function proforma_results(p::REoptInputs, d::Dict)
     end
 
     # calculate Storage o+m costs, incentives, and depreciation
-    if "ElectricStorage" in keys(d) && d["ElectricStorage"]["size_kw"] > 0
-        # TODO handle other types of storage
-        storage = p.s.storage.attr["ElectricStorage"]
-        total_kw = d["ElectricStorage"]["size_kw"]
-        total_kwh = d["ElectricStorage"]["size_kwh"]
-        capital_cost = total_kw * storage.installed_cost_per_kw + total_kwh * storage.installed_cost_per_kwh
-        battery_replacement_year = storage.battery_replacement_year
-        battery_replacement_cost = -1 * ((total_kw * storage.replace_cost_per_kw) + (
-                    total_kwh * storage.replace_cost_per_kwh))
-        m.om_series += [yr != battery_replacement_year ? 0 : battery_replacement_cost for yr in 1:years]
-
-        # storage only has cbi in the API
-        cbi = total_kw * storage.total_rebate_per_kw + total_kwh * storage.total_rebate_per_kwh
-        m.total_ibi_and_cbi += cbi
-
-        # ITC
-        federal_itc_basis = capital_cost  # bug in v1 subtracted cbi from capital_cost here
-        federal_itc_amount = storage.total_itc_fraction * federal_itc_basis
-        m.federal_itc += federal_itc_amount
-
-        # Depreciation
-        if storage.macrs_option_years in [5 ,7]
-            depreciation_schedule = get_depreciation_schedule(p, storage, federal_itc_basis)
-            m.total_depreciation += depreciation_schedule
+    if "ElectricStorage" in keys(d)
+        if typeof(d["ElectricStorage"]) <: AbstractArray && length(d["ElectricStorage"]) == 1
+            d["ElectricStorage"] = d["ElectricStorage"][1]
+        elseif !(typeof(d["ElectricStorage"]) <: AbstractDict)
+            throw(@error("Pro forma results for a REopt solution with multiple ElectricStorage is not yet supported"))
+        end
+        elec_stor_name = get(d["ElectricStorage"],"name","ElectricStorage")
+        if d["ElectricStorage"]["size_kw"] > 0
+            # TODO handle other types of storage
+            storage = p.s.storage.attr[elec_stor_name]
+            total_kw = d["ElectricStorage"]["size_kw"]
+            total_kwh = d["ElectricStorage"]["size_kwh"]
+            capital_cost = total_kw * storage.installed_cost_per_kw + total_kwh * storage.installed_cost_per_kwh
+            battery_replacement_year = storage.battery_replacement_year
+            battery_replacement_cost = -1 * ((total_kw * storage.replace_cost_per_kw) + (
+                        total_kwh * storage.replace_cost_per_kwh))
+            m.om_series += [yr != battery_replacement_year ? 0 : battery_replacement_cost for yr in 1:years]
+
+            # storage only has cbi in the API
+            cbi = total_kw * storage.total_rebate_per_kw + total_kwh * storage.total_rebate_per_kwh
+            m.total_ibi_and_cbi += cbi
+
+            # ITC
+            federal_itc_basis = capital_cost  # bug in v1 subtracted cbi from capital_cost here
+            federal_itc_amount = storage.total_itc_fraction * federal_itc_basis
+            m.federal_itc += federal_itc_amount
+
+            # Depreciation
+            if storage.macrs_option_years in [5, 7]
+                depreciation_schedule = get_depreciation_schedule(p, storage, federal_itc_basis)
+                m.total_depreciation += depreciation_schedule
+            end
         end
     end
 
diff --git a/src/results/pv.jl b/src/results/pv.jl
index 4c5774faa..02cda39ac 100644
--- a/src/results/pv.jl
+++ b/src/results/pv.jl
@@ -66,6 +66,8 @@ function add_pv_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 		PVPerUnitSizeOMCosts = p.third_party_factor * p.om_cost_per_kw[t] * p.pwf_om * m[Symbol("dvSize"*_n)][t]
 		r["lifecycle_om_cost_after_tax"] = round(value(PVPerUnitSizeOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
         r["lcoe_per_kwh"] = calculate_lcoe(p, r, get_pv_by_name(t, p.s.pvs))
+        r["year_one_om_cost_before_tax"] = round(value(PVPerUnitSizeOMCosts) / p.pwf_om, digits=0)
+        r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)][t]) * get_pv_by_name(t, p.s.pvs).installed_cost_per_kw, digits=3)
         d[t] = r
 	end
     nothing
diff --git a/src/results/steam_turbine.jl b/src/results/steam_turbine.jl
index e25aa37a2..48fe2b40b 100644
--- a/src/results/steam_turbine.jl
+++ b/src/results/steam_turbine.jl
@@ -25,6 +25,7 @@ function add_steam_turbine_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dic
     r = Dict{String, Any}()
 
 	r["size_kw"] = round(value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.steam_turbine)), digits=3)
+	r["initial_capital_cost"] = round(value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.steam_turbine)) * p.s.steam_turbine.installed_cost_per_kw, digits=3)
     @expression(m, Year1SteamTurbineThermalConsumptionKWH,
 		p.hours_per_time_step * sum(m[Symbol("dvThermalToSteamTurbine"*_n)][tst,q,ts] for tst in p.techs.can_supply_steam_turbine, q in p.heating_loads, ts in p.time_steps))
     r["annual_thermal_consumption_mmbtu"] = round(value(Year1SteamTurbineThermalConsumptionKWH) / KWH_PER_MMBTU, digits=5)
@@ -50,7 +51,7 @@ function add_steam_turbine_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dic
 	r["electric_to_grid_series_kw"] = round.(value.(SteamTurbinetoGrid), digits=3)
 	if !isempty(p.s.storage.types.elec)
 		@expression(m, SteamTurbinetoBatt[ts in p.time_steps],
-			sum(m[Symbol("dvProductionToStorage"*_n)]["ElectricStorage",t,ts] for t in p.techs.steam_turbine))
+			sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.steam_turbine, b in p.s.storage.types.elec))
 	else
 		SteamTurbinetoBatt = zeros(length(p.time_steps))
 	end
@@ -99,7 +100,6 @@ function add_steam_turbine_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dic
         @expression(m, SteamTurbineToProcessHeatKW[ts in p.time_steps], 0.0)
     end
     r["thermal_to_process_heat_load_series_mmbtu_per_hour"] = round.(value.(SteamTurbineToProcessHeatKW ./ KWH_PER_MMBTU), digits=5)
-
 	
 	d["SteamTurbine"] = r
 	nothing
diff --git a/src/results/thermal_storage.jl b/src/results/thermal_storage.jl
index 48e9f607c..602269cb4 100644
--- a/src/results/thermal_storage.jl
+++ b/src/results/thermal_storage.jl
@@ -20,6 +20,7 @@ function add_hot_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict,
     r = Dict{String, Any}()
     size_kwh = round(value(m[Symbol("dvStorageEnergy"*_n)][b]), digits=3)
     r["size_gal"] = round(size_kwh / kwh_per_gal, digits=0)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvStorageEnergy"*_n)][b]) * p.s.storage.attr[b].installed_cost_per_kwh, digits=3)
 
     if size_kwh != 0
     	soc = (m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps)
@@ -93,12 +94,13 @@ function add_cold_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict
     Note: the node number is an empty string if evaluating a single `Site`.
     =#
 
-    kwh_per_gal = get_kwh_per_gal(p.s.storage.attr["ColdThermalStorage"].hot_water_temp_degF,
-                                    p.s.storage.attr["ColdThermalStorage"].cool_water_temp_degF)
+    kwh_per_gal = get_kwh_per_gal(p.s.storage.attr[b].hot_water_temp_degF,
+                                    p.s.storage.attr[b].cool_water_temp_degF)
     
     r = Dict{String, Any}()
     size_kwh = round(value(m[Symbol("dvStorageEnergy"*_n)][b]), digits=3)
     r["size_gal"] = round(size_kwh / kwh_per_gal, digits=0)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvStorageEnergy"*_n)][b]) * p.s.storage.attr[b].installed_cost_per_kwh, digits=3)
 
     if size_kwh != 0
     	soc = (m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps)
diff --git a/src/results/wind.jl b/src/results/wind.jl
index 0e89e19dc..7c549d1de 100644
--- a/src/results/wind.jl
+++ b/src/results/wind.jl
@@ -26,7 +26,7 @@ function add_wind_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 	r["production_factor_series"] = Vector(p.production_factor[t, :])
 	per_unit_size_om = @expression(m, p.third_party_factor * p.pwf_om * m[:dvSize][t] * p.om_cost_per_kw[t])
 
-	r["size_kw"] = round(value(m[:dvSize][t]), digits=2)
+	r["size_kw"] = round(value(m[Symbol("dvSize"*_n)][t]), digits=2)
 	r["lifecycle_om_cost_after_tax"] = round(value(per_unit_size_om) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
 	r["year_one_om_cost_before_tax"] = round(value(per_unit_size_om) / (p.pwf_om * p.third_party_factor), digits=0)
 
@@ -64,6 +64,7 @@ function add_wind_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 	r["annual_energy_produced_kwh"] = round(value(AvgWindProd), digits=0)
 
     r["lcoe_per_kwh"] = calculate_lcoe(p, r, p.s.wind)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)][t]) * p.s.wind.installed_cost_per_kw, digits=3)
 	d[t] = r
     nothing
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0d355fae0..fddd28a12 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,6 +23,21 @@ elseif "CPLEX" in ARGS
     end
 else  # run HiGHS tests
     @testset verbose=true "REopt test set using HiGHS solver" begin
+        @testset "Multiple Electric Storage" begin
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            r = run_reopt(model, "./scenarios/multiple_elec_storage.json")
+            @test r["PV"]["size_kw"] ≈ 216.6667 atol=0.01
+            @test r["Financial"]["lcc"] ≈ 1.2391786e7 rtol=1e-5
+            for stor in r["ElectricStorage"]
+                if stor["name"] == "CheapElecStor"
+                    @test stor["size_kw"] ≈ 49.0 atol=0.1
+                    @test stor["size_kwh"] ≈ 83.3 atol=0.1
+                else
+                    @test stor["size_kw"] ≈ 0 atol=0.1
+                    @test stor["size_kwh"] ≈ 0 atol=0.1
+                end
+            end
+        end
         @testset "Prevent simultaneous charge and discharge" begin
             logger = SimpleLogger()
             results = nothing
@@ -656,6 +671,66 @@ else  # run HiGHS tests
             end
         end
 
+        @testset "Electric Storage O&M" begin
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            data = JSON.parsefile("./scenarios/storage_om.json")
+
+            data["ElectricStorage"]["om_cost_per_kw"] = 10
+            data["ElectricStorage"]["om_cost_per_kwh"] = 0
+
+            inputs = REoptInputs(data)
+            results1 = run_reopt(model, inputs)
+
+            @test results1["ElectricStorage"]["year_one_om_cost_before_tax"] ≈ 400 atol=1
+            @test results1["ElectricStorage"]["lifecycle_om_cost_after_tax"] ≈ 6272 atol=1 
+
+            data["ElectricStorage"]["om_cost_per_kw"] = 10
+            data["ElectricStorage"]["om_cost_per_kwh"] = 5
+
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            inputs = REoptInputs(data)
+            results2 = run_reopt(model, inputs)
+
+            @test results2["ElectricStorage"]["year_one_om_cost_before_tax"] ≈ 800 atol=1
+            @test results2["ElectricStorage"]["lifecycle_om_cost_after_tax"] ≈ 12543 atol=1 
+            @test results2["Financial"]["lcc"] ≈ results1["Financial"]["lcc"] + 12543 - 6272 atol=1
+        end
+
+        @testset "Electric Storage Self-Discharge" begin
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            data = JSON.parsefile("./scenarios/storage_om.json")
+
+            data["ElectricStorage"]["soc_init_fraction"] = 1.0
+            data["ElectricStorage"]["soc_min_fraction"] = 0.0
+            data["ElectricStorage"]["soc_based_per_ts_self_discharge_fraction"] = 0.0025/24
+
+            s = Scenario(data)
+            inputs = REoptInputs(s)
+            results_soc_based_self_discharge = run_reopt(model, inputs)
+
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            data["ElectricStorage"]["soc_based_per_ts_self_discharge_fraction"] = 0.0
+
+            s = Scenario(data)
+            inputs = REoptInputs(s)
+            results_no_self_discharge = run_reopt(model, inputs)
+
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            data["ElectricStorage"]["capacity_based_per_ts_self_discharge_fraction"] = 0.0025/24
+
+            s = Scenario(data)
+            inputs = REoptInputs(s)
+            results_capacity_based_self_discharge = run_reopt(model, inputs)
+
+            @test results_soc_based_self_discharge["ElectricStorage"]["year_one_om_cost_before_tax"] ≈ 0 atol=0.01
+            @test results_soc_based_self_discharge["ElectricStorage"]["lifecycle_om_cost_after_tax"] ≈ 0 atol=0.01 
+            @test sum(results_no_self_discharge["ElectricStorage"]["storage_to_load_series_kw"]) - 
+                sum(results_soc_based_self_discharge["ElectricStorage"]["storage_to_load_series_kw"]) ≈ 15.382 atol=0.01 
+
+            calculated_output = round((80 - (80 * 0.0025/24 * 8760)) * sqrt(0.975) * 0.96, digits = 3)
+            @test sum(results_capacity_based_self_discharge["ElectricStorage"]["storage_to_load_series_kw"]) ≈ calculated_output atol=0.01 
+        end
+        
         @testset "Heating loads and addressable load fraction" begin
             # Default LargeOffice CRB with SpaceHeatingLoad and DomesticHotWaterLoad are served by ExistingBoiler
             m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
@@ -2925,4 +3000,4 @@ else  # run HiGHS tests
 
         end
     end
-end
\ No newline at end of file
+end 
\ No newline at end of file
diff --git a/test/scenarios/ghp_inputs.json b/test/scenarios/ghp_inputs.json
index 43b816527..abd01faad 100644
--- a/test/scenarios/ghp_inputs.json
+++ b/test/scenarios/ghp_inputs.json
@@ -59,11 +59,11 @@
     "ColdThermalStorage": {
         "min_gal": 10,
         "max_gal": 10,
-        "thermal_decay_rate_fraction": 0.0
+        "capacity_based_per_ts_thermal_decay_fraction": 0.0
     },
     "HotThermalStorage": {
         "min_gal": 10,
         "max_gal": 10,
-        "thermal_decay_rate_fraction": 0.0
+        "capacity_based_per_ts_thermal_decay_fraction": 0.0
     }
 }
\ No newline at end of file
diff --git a/test/scenarios/heat_cool_energy_balance_inputs.json b/test/scenarios/heat_cool_energy_balance_inputs.json
index 6d8af917e..30aecc394 100644
--- a/test/scenarios/heat_cool_energy_balance_inputs.json
+++ b/test/scenarios/heat_cool_energy_balance_inputs.json
@@ -74,13 +74,13 @@
         "min_gal": 20000.0,
         "max_gal": 20000.0,
         "internal_efficiency_fraction": 0.97,
-        "thermal_decay_rate_fraction": 0.004
+        "capacity_based_per_ts_thermal_decay_fraction": 0.004
     },
     "ColdThermalStorage":{
         "min_gal": 30000.0,
         "max_gal": 30000.0,
         "internal_efficiency_fraction": 0.97,
-        "thermal_decay_rate_fraction": 0.004,
+        "capacity_based_per_ts_thermal_decay_fraction": 0.004,
         "soc_init_fraction": 0.1
     },
     "AbsorptionChiller": {
diff --git a/test/scenarios/multiple_elec_storage.json b/test/scenarios/multiple_elec_storage.json
new file mode 100644
index 000000000..e0704c4e1
--- /dev/null
+++ b/test/scenarios/multiple_elec_storage.json
@@ -0,0 +1,68 @@
+{
+    "Site": {
+        "longitude": -118.1164613,
+        "latitude": 34.5794343,
+        "roof_squarefeet": 5000.0,
+        "land_acres": 1.0,
+        "node": 3
+    },
+    "ElectricUtility" : {
+        "co2_from_avert" : true
+    },
+    "PV": {
+        "macrs_bonus_fraction": 0.4,
+        "installed_cost_per_kw": 2000.0,
+        "tilt": 34.579,
+        "degradation_fraction": 0.005,
+        "macrs_option_years": 5,
+        "federal_itc_fraction": 0.3,
+        "module_type": 0,
+        "array_type": 1,
+        "om_cost_per_kw": 16.0,
+        "macrs_itc_reduction": 0.5,
+        "azimuth": 180.0,
+        "federal_rebate_per_kw": 350.0
+    },
+    "ElectricLoad": {
+        "doe_reference_name": "RetailStore",
+        "annual_kwh": 10000000.0,
+        "year": 2017
+    },
+    "ElectricStorage": [
+        {
+            "name": "CheapElecStor",
+            "total_rebate_per_kw": 100.0,
+            "macrs_option_years": 5,
+            "can_grid_charge": true,
+            "macrs_bonus_fraction": 0.4,
+            "replace_cost_per_kw": 460.0,
+            "replace_cost_per_kwh": 230.0,
+            "installed_cost_per_kw": 1000.0,
+            "installed_cost_per_kwh": 500.0,
+            "total_itc_fraction": 0.0
+        },
+        {
+            "total_rebate_per_kw": 100.0,
+            "macrs_option_years": 5,
+            "can_grid_charge": true,
+            "macrs_bonus_fraction": 0.4,
+            "replace_cost_per_kw": 460.0,
+            "replace_cost_per_kwh": 230.0,
+            "installed_cost_per_kw": 1500.0,
+            "installed_cost_per_kwh": 500.0,
+            "total_itc_fraction": 0.0
+        }
+    ],
+    "ElectricTariff": {
+        "urdb_label": "5ed6c1a15457a3367add15ae"
+    },
+    "Financial": {
+        "elec_cost_escalation_rate_fraction": 0.026,
+        "offtaker_discount_rate_fraction": 0.081,
+        "owner_discount_rate_fraction": 0.081,
+        "analysis_years": 20,
+        "offtaker_tax_rate_fraction": 0.4,
+        "owner_tax_rate_fraction": 0.4,
+        "om_cost_escalation_rate_fraction": 0.025
+    }
+}
diff --git a/test/scenarios/re_emissions_with_thermal.json b/test/scenarios/re_emissions_with_thermal.json
index fbc0346f1..747f1792e 100644
--- a/test/scenarios/re_emissions_with_thermal.json
+++ b/test/scenarios/re_emissions_with_thermal.json
@@ -79,7 +79,7 @@
     "soc_min_fraction": 0.1 ,
     "soc_init_fraction": 0.5 ,
     "installed_cost_per_gal": 3.0 ,
-    "thermal_decay_rate_fraction": 0.004 ,
+    "capacity_based_per_ts_thermal_decay_fraction": 0.004 ,
     "om_cost_per_gal": 0.0 ,
     "macrs_option_years": 0 ,
     "macrs_bonus_fraction": 0.0
diff --git a/test/scenarios/storage_om.json b/test/scenarios/storage_om.json
new file mode 100644
index 000000000..c188b51b8
--- /dev/null
+++ b/test/scenarios/storage_om.json
@@ -0,0 +1,46 @@
+{
+    "Site": {
+        "longitude": -118.1164613,
+        "latitude": 34.5794343,
+        "roof_squarefeet": 5000.0,
+        "land_acres": 1.0
+    },
+    "ElectricLoad": {
+        "doe_reference_name": "RetailStore",
+        "annual_kwh": 10000000.0,
+        "year": 2017
+    },
+    "ElectricStorage": {
+        "min_kw": 40,
+        "max_kw": 40,
+        "min_kwh": 80,
+        "max_kwh": 80,   
+        "can_grid_charge": false,
+        "replace_cost_per_kw": 460.0,
+        "replace_cost_per_kwh": 230.0,
+        "installed_cost_per_kw": 1000.0,
+        "installed_cost_per_kwh": 500.0,
+        "total_itc_fraction": 0.0,
+        "macrs_option_years": 0,
+        "macrs_bonus_fraction": 0,
+        "soc_init_fraction": 1.0,
+        "soc_min_fraction": 0.0,
+        "can_export_to_grid": true
+    },
+    "ElectricTariff": {
+        "urdb_label": "5ed6c1a15457a3367add15ae"
+    },
+    "PV": {
+        "min_kw": 0.0,
+        "max_kw": 0.0
+    },
+    "Financial": {
+        "elec_cost_escalation_rate_fraction": 0.026,
+        "offtaker_discount_rate_fraction": 0.05,
+        "owner_discount_rate_fraction": 0.05,
+        "analysis_years": 20,
+        "offtaker_tax_rate_fraction": 0,
+        "owner_tax_rate_fraction": 0,
+        "om_cost_escalation_rate_fraction": 0.025
+    }
+}