diff --git a/src/model/resources/thermal_storage/thermal_storage.jl b/src/model/resources/thermal_storage/thermal_storage.jl
index 8983e63122..7a101b3fc5 100644
--- a/src/model/resources/thermal_storage/thermal_storage.jl
+++ b/src/model/resources/thermal_storage/thermal_storage.jl
@@ -109,7 +109,9 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 	dfTS = inputs["dfTS"]
 	RH = get_resistive_heating(inputs)
 
+
 	by_rid(rid, sym) = by_rid_df(rid, sym, dfTS)
+	load_fuel_data!(inputs, setup)
 
 	@variables(EP, begin
 	# Thermal core variables
@@ -136,7 +138,7 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 
 	# Variable cost of core operation
 	# Variable cost at timestep t for thermal core y
-	@expression(EP, eCVar_Core[y in TS, t=1:T], inputs["omega"][t] * by_rid(y, :Var_OM_Cost_per_MWh_th) * vCP[y,t])
+	@expression(EP, eCVar_Core[y in TS, t=1:T], inputs["omega"][t] * (by_rid(y, :Var_OM_Cost_per_MWh_th) + inputs["TS_C_Fuel_per_MWh"][y][t]) * vCP[y,t])
 	# Variable cost from all thermal cores at timestep t)
 	@expression(EP, eTotalCVarCoreT[t=1:T], sum(eCVar_Core[y,t] for y in TS))
 	# Total variable cost for all thermal cores
@@ -295,9 +297,66 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 		EP[:eCapResMarBalance] += eCapResMarBalanceFusionAdjustment
 	end
 
+	# add emissions
+	thermal_core_emissions!(EP, inputs, setup)
+
 	return EP
 end
 
+function load_fuel_data!(inputs::Dict, setup::Dict)
+
+	dfTS = inputs["dfTS"]
+	TS = inputs["TS"]
+	NONFUS = get_nonfus(inputs)
+	N = nrow(dfTS)
+
+	# for unit commitment decisions
+	if setup["UCommit"]>=1
+		# Convert to $ million/GW with objective function in millions
+		if setup["ParameterScale"] ==1	
+			dfTS[!,:Start_Cost_per_MW] = dfTS[!,:Start_Cost_per_MW]/ModelScalingFactor 
+		end
+
+		# Fuel consumed on start-up (million BTUs per MW per start) if unit commitment is modelled
+		start_fuel = convert(Array{Float64}, collect(skipmissing(dfTS[!,:Start_Fuel_MMBTU_per_MW])))
+		# Fixed cost per start-up ($ per MW per start) if unit commitment is modelled
+		start_cost = convert(Array{Float64}, collect(skipmissing(dfTS[!,:Start_Cost_per_MW])))
+		inputs["TS_C_Start"] = Dict()
+		dfTS[!,:CO2_per_Start] = zeros(Float64, N)
+	end
+
+	# Heat rate of all resources (million BTUs/MWh)
+	heat_rate = convert(Array{Float64}, collect(skipmissing(dfTS[!,:Heat_Rate_MMBTU_per_MWh])))
+	# Fuel used by each resource
+	fuel_type = collect(skipmissing(dfTS[!,:Fuel]))
+	# Maximum fuel cost in $ per MWh and CO2 emissions in tons per MWh
+	inputs["TS_C_Fuel_per_MWh"] = Dict()
+	dfTS[!,:CO2_per_MWh] = zeros(Float64,N)
+
+	for g in 1:N 
+		#calculate fuel costs
+		inputs["TS_C_Fuel_per_MWh"][dfTS[g, :R_ID]] = inputs["fuel_costs"][fuel_type[g]] .* heat_rate[g]
+		#calculate fuel emissions
+		dfTS[g, :CO2_per_MWh] = inputs["fuel_CO2"][fuel_type[g]] .* heat_rate[g]
+		if setup["ParameterScale"] ==1
+			dfTS[g,:CO2_per_MWh] = dfTS[g,:CO2_per_MWh] * ModelScalingFactor
+		end
+
+		# add start up costs and emissions for committed thermal cores.
+		if dfTS[g, :R_ID] in intersect(inputs["THERM_COMMIT"])
+			inputs["TS_C_Start"][dfTS[g, :R_ID]] = by_rid_df(g, :Cap_Size, dfTS) .* (inputs["fuel_costs"][fuel_type[g]] .* start_fuel[g] .+ start_cost[g])
+
+			dfTS[g, :CO2_per_Start] = by_rid_df(g, :Cap_Size, dfTS) * (inputs["fuel_CO2"][fuel_type[g]] * start_fuel[g])
+
+			#scale appropriately
+			if setup["ParameterScale"] == 1
+				dfTS[g, :CO2_per_Start] = dfTS[g, :CO2_per_Start] * ModelScalingFactor
+			end
+		end
+	end
+end
+
+
 function fusion_max_cap_constraint!(EP::Model, inputs::Dict, setup::Dict)
 
 	dfGen = inputs["dfGen"]
@@ -485,7 +544,7 @@ function thermal_core_constraints!(EP::Model, inputs::Dict, setup::Dict)
 	INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]
 
 	COMMIT = intersect(inputs["THERM_COMMIT"], NONFUS)
-	NON_COMMIT = intersect(inputs["THERM_NO_COMMIT"])
+	NON_COMMIT = intersect(inputs["THERM_NO_COMMIT"], NONFUS)
 
 	# constraints for generators not subject to UC
 	if !isempty(NON_COMMIT)
@@ -521,7 +580,11 @@ function thermal_core_constraints!(EP::Model, inputs::Dict, setup::Dict)
 		# shutdown event variable
 		@variable(EP, vCSHUT[y in COMMIT, t=1:T] >= 0)
 
-		### TODO: STARTUP COSTS ???????? ###
+		### Add startup costs ###
+		@expression(EP, eCStartTS[y in COMMIT, t=1:T], (inputs["omega"][t] * inputs["TS_C_Start"][y][t] * vCSTART[y, t]))
+		@expression(EP, eTotalCStartTST[t=1:T], sum(eCStartTS[y,t] for y in COMMIT))
+		@expression(EP, eTotalCStartTS, sum(eTotalCStartTST[t] for t=1:T))
+		EP[:eObj] += eTotalCStartTS
 
 		## Declaration of integer/binary variables
 		if setup["UCommit"] == 1 # Integer UC constraints
@@ -655,6 +718,39 @@ function maintenance_constraints!(EP::Model, inputs::Dict, setup::Dict)
 		(EP[:vCCAP][y] - by_rid(y,:Cap_Size) * EP[:vFMDOWN][y,t]) * dfGen[y,:Eff_Down] * by_rid(y,:Recirc_Pass))
 end
 
+function thermal_core_emissions!(EP::Model, inputs::Dict, setup::Dict)
+
+	dfTS = inputs["dfTS"]
+	dfGen = inputs["dfGen"]
+
+	TS = inputs["TS"]	# R_IDs of resources with thermal storage
+	G = inputs["G"]		# R_IDs of all resources
+	T = inputs["T"]     # Number of time steps (hours)
+	Z = inputs["Z"]     # Number of zones
+	FUS = get_fus(inputs) #FUS generators
+	NONFUS = get_nonfus(inputs)	#NONFUS generators
+	by_rid(rid, sym) = by_rid_df(rid, sym, dfTS)
+	
+	@expression(EP, eEmissionsByPlantTS[y = 1:G, t = 1:T],
+
+		if y ∉ TS
+			0
+		elseif y in FUS
+			by_rid(y, :CO2_per_MWh) * EP[:vCP][y, t] + by_rid(y, :CO2_per_Start) * EP[:vFSTART][y, t]
+		elseif y in intersect(NONFUS, inputs["THERM_COMMIT"])
+			by_rid(y, :CO2_per_MWh) * EP[:vCP][y, t] + by_rid(y, :CO2_per_Start) * EP[:vCSTART][y, t]
+		else
+			by_rid(y, :CO2_per_MWh) * EP[:vCP][y,t]
+		end
+	)
+
+	@expression(EP, eEmissionsByZoneTS[z=1:Z, t=1:T], sum(eEmissionsByPlantTS[y,t] for y in intersect(TS, dfGen[(dfGen[!,:Zone].==z),:R_ID])))
+
+	EP[:eEmissionsByPlant] += eEmissionsByPlantTS
+	EP[:eEmissionsByZone] += eEmissionsByZoneTS
+
+end
+
 
 function sanity_check_maintenance(MAINTENANCE::Vector{Int}, setup::Dict)
 	ow = setup["OperationWrapping"]
diff --git a/src/write_outputs/write_thermal_storage.jl b/src/write_outputs/write_thermal_storage.jl
index 484b2495db..37f6a9600a 100644
--- a/src/write_outputs/write_thermal_storage.jl
+++ b/src/write_outputs/write_thermal_storage.jl
@@ -15,8 +15,55 @@ received this license file.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 
+function write_core_capacities(EP::Model, inputs::Dict, filename::AbstractString, msf)
 
-function write_core_behaviors(EP::Model, inputs::Dict, symbol::Symbol, SET::Vector{Int}, filename::AbstractString)
+	# Capacity decisions
+	dfGen = inputs["dfGen"]
+	dfTS = inputs["dfTS"]
+	T = inputs["T"]
+
+	# load capacity power
+	TSResources = dfTS[!,:Resource]
+	TSG = length(TSResources)
+	corecappower = zeros(TSG)
+	for i in 1:TSG
+		corecappower[i] = first(value.(EP[:vCCAP][dfTS[i,:R_ID]]))
+	end
+
+	# load capacity energy
+	corecapenergy = zeros(TSG)
+	for i in 1:TSG
+		corecapenergy[i] = first(value.(EP[:vTSCAP][dfTS[i,:R_ID]]))
+	end
+
+	# load rh capacity
+	rhcapacity = zeros(TSG)
+	for i in 1:TSG
+		if dfTS[i, :RH] == 1
+			rhcapacity[i] = first(value.(EP[:vRHCAP][dfTS[i,:R_ID]]))
+		end
+	end
+
+	# create data frame
+	dfCoreCap = DataFrame(
+		Resource = TSResources, 
+		Zone = dfTS[!,:Zone],
+		CorePowerCap = corecappower[:],
+		TSEnergyCap = corecapenergy[:],
+		RHPowerCap = rhcapacity[:]
+	)
+
+	# adjust files and write
+	dfCoreCap.CorePowerCap = dfCoreCap.CorePowerCap * msf
+	dfCoreCap.TSEnergyCap = dfCoreCap.TSEnergyCap * msf
+	dfCoreCap.RHPowerCap = dfCoreCap.RHPowerCap * msf
+	CSV.write(filename, dfCoreCap)
+
+	return dfCoreCap
+
+end
+
+function write_core_commitments(EP::Model, inputs::Dict, SET::Vector{Int},symbol::Symbol, filename::AbstractString)
 	dfTS = inputs["dfTS"]
 	T = inputs["T"]
 
@@ -39,7 +86,7 @@ function write_core_behaviors(EP::Model, inputs::Dict, symbol::Symbol, SET::Vect
 	return df
 end
 
-function write_scaled_values(EP::Model, inputs::Dict, symbol::Symbol, SET::Vector{Int}, filename::AbstractString, msf)
+function write_scaled_values(EP::Model, inputs::Dict, SET::Vector{Int}, symbol::Symbol, filename::AbstractString, msf)
 	dfTS = inputs["dfTS"]
 	T = inputs["T"]
 
@@ -62,106 +109,100 @@ function write_scaled_values(EP::Model, inputs::Dict, symbol::Symbol, SET::Vecto
 	return df
 end
 
-function write_thermal_storage_system_max_dual(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_thermal_storage_system_max_dual(EP::Model, inputs::Dict, setup::Dict, filename::AbstractString, msf)
 	dfTS = inputs["dfTS"]
+	FUS = dfTS[dfTS.FUS .== 1, :R_ID]
 
-	FIRST_ROW = 1
-	if dfTS[FIRST_ROW, :System_Max_Cap_MWe_net] >= 0
-		val = -1*dual.(EP[:cCSystemTot])
-		val *= setup["ParameterScale"] == 1 ? ModelScalingFactor : 1.0
-		df = DataFrame(:System_Max_Cap_MW_th_dual => val)
-		filename = joinpath(path, "System_Max_TS_Cap_dual.csv")
+	if !isempty(FUS)
+		FIRST_ROW = 1
+		if dfTS[FIRST_ROW, :System_Max_Cap_MWe_net] >= 0
+			val = -1*dual.(EP[:cCSystemTot])
+			val *= msf
+			df = DataFrame(:System_Max_Cap_MW_th_dual => val)
+			CSV.write(filename, dftranspose(df, false), writeheader=false)
+		end
+	end	
+end
+
+function write_thermal_storage_capacity_duals(EP::Model, inputs::Dict, setup::Dict, filename::AbstractString, msf)
+	dfTS = inputs["dfTS"]
+	NONFUS = dfTS[dfTS.FUS .== 0, :R_ID]
+
+	if !isempty(NONFUS)
+		HAS_MAX_LIMIT = dfTS[by_rid_df(NONFUS, :Max_Core_Power_Capacity, dfTS) .> 0, :R_ID]
+		resources = by_rid_df(HAS_MAX_LIMIT, :Resource, dfTS)
+		n_max = length(HAS_MAX_LIMIT)
+		vals = zeros(n_max)
+		for i in 1:n_max
+			vals[i] = -1 * dual.(EP[:cCoreMaxCapacity][i]) * msf
+		end
+		df = DataFrame(
+			Resource = resources,
+			R_ID = HAS_MAX_LIMIT,
+			Dual = vals
+		)
 		CSV.write(filename, dftranspose(df, false), writeheader=false)
 	end
 end
 
-
 @doc raw"""
 	write_capacity(path::AbstractString, inputs::Dict, setup::Dict, EP::Model))
 
 Function for writing the diferent capacities for the different generation technologies (starting capacities or, existing capacities, retired capacities, and new-built capacities).
 """
 function write_thermal_storage(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	# Capacity decisions
+
+	### LOAD DATASETS AND PREPARE SCALING FACTOR
 	dfGen = inputs["dfGen"]
 	dfTS = inputs["dfTS"]
 	T = inputs["T"]
-
-
-	TSResources = dfTS[!,:Resource]
-	TSG = length(TSResources)
-	corecappower = zeros(TSG)
-	for i in 1:TSG
-		corecappower[i] = first(value.(EP[:vCCAP][dfTS[i,:R_ID]]))
-	end
-
-	corecapenergy = zeros(TSG)
-	for i in 1:TSG
-		corecapenergy[i] = first(value.(EP[:vTSCAP][dfTS[i,:R_ID]]))
-	end
-
-	rhcapacity = zeros(TSG)
-	for i in 1:TSG
-		if dfTS[i, :RH] == 1
-			rhcapacity[i] = first(value.(EP[:vRHCAP][dfTS[i,:R_ID]]))
-		else
-			rhcapacity[i] = 0
-		end
-	end
-
-	dfCoreCap = DataFrame(
-		Resource = TSResources, Zone = dfTS[!,:Zone],
-		CorePowerCap = corecappower[:],
-		TSEnergyCap = corecapenergy[:],
-		RHPowerCap = rhcapacity[:]
-	)
-
-	# set a single scalar to avoid future branching
 	msf = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-	dfCoreCap.CorePowerCap = dfCoreCap.CorePowerCap * msf
-	dfCoreCap.TSEnergyCap = dfCoreCap.TSEnergyCap * msf
-	dfCoreCap.RHPowerCap = dfCoreCap.RHPowerCap * msf
-	CSV.write(joinpath(path,"TS_capacity.csv"), dfCoreCap)
 
+	### WRITE CORE CAPACITY DECISIONS ###
+	dfCoreCap = write_core_capacities(EP, inputs, joinpath(path,"TS_capacity.csv"), msf)
+
+	### LOAD RELEVANT SETS ###
 	THERMAL_STORAGE = dfTS.R_ID
 	RH = dfTS[dfTS.RH .==1, :R_ID]
 	FUS = dfTS[dfTS.FUS .== 1, :R_ID]
 	NONFUS = dfTS[dfTS.FUS .== 0, :R_ID]
+
 	### CORE POWER TIME SERIES ###
-	dfCorePwr = write_scaled_values(EP, inputs, :vCP, THERMAL_STORAGE, joinpath(path, "TS_CorePwr.csv"), msf)
+	dfCorePwr = write_scaled_values(EP, inputs, THERMAL_STORAGE, :vCP, joinpath(path, "TS_CorePwr.csv"), msf)
 
 	### THERMAL SOC TIME SERIES ###
-	dfTSOC = write_scaled_values(EP, inputs, :vTS, THERMAL_STORAGE, joinpath(path, "TS_SOC.csv"), msf)
+	dfTSOC = write_scaled_values(EP, inputs, THERMAL_STORAGE, :vTS, joinpath(path, "TS_SOC.csv"), msf)
 
 	### RESISTIVE HEATING TIME SERIES ### 
 	if !isempty(RH)
-		dfRH = write_scaled_values(EP, inputs, :vRH, RH, joinpath(path, "TS_RH.csv"), msf)
+		dfRH = write_scaled_values(EP, inputs, RH, :vRH, joinpath(path, "TS_RH.csv"), msf)
 	end
 
 	### FUSION SPECIFIC OUTPUTS ###
 	if !isempty(FUS)
 		### RECIRCULATING POWER TIME SERIES ###
-		dfRecirc = write_scaled_values(EP, inputs, :eTotalRecircFus, FUS, joinpath(path, "TS_Recirc.csv"), msf)
+		dfRecirc = write_scaled_values(EP, inputs, FUS, :eTotalRecircFus, joinpath(path, "TS_Recirc.csv"), msf)
 
 		### CORE STARTS, SHUTS, COMMITS, and MAINTENANCE TIMESERIES ###
-		dfFStart = write_core_behaviors(EP, inputs, :vFSTART, FUS, joinpath(path, "TS_FUS_start.csv"))
-		dfFShut = write_core_behaviors(EP, inputs, :vFSHUT, FUS, joinpath(path, "TS_FUS_shut.csv"))
-		dfFCommit = write_core_behaviors(EP, inputs, :vFCOMMIT, FUS, joinpath(path, "TS_FUS_commit.csv"))
+		dfFStart = write_core_commitments(EP, inputs, FUS, :vFSTART, joinpath(path, "TS_FUS_start.csv"))
+		dfFShut = write_core_commitments(EP, inputs, FUS, :vFSHUT, joinpath(path, "TS_FUS_shut.csv"))
+		dfFCommit = write_core_commitments(EP, inputs, FUS, :vFCOMMIT, joinpath(path, "TS_FUS_commit.csv"))
 
 		if setup["OperationWrapping"] == 0 && !isempty(get_maintenance(inputs))
-			dfMaint = write_core_behaviors(EP, inputs, :vFMDOWN, FUS,  joinpath(path, "TS_FUS_maint.csv"))
-			dfMShut = write_core_behaviors(EP, inputs, :vFMSHUT, FUS, joinpath(path, "TS_FUS_maintshut.csv"))
+			dfMaint = write_core_commitments(EP, inputs, FUS, :vFMDOWN, joinpath(path, "TS_FUS_maint.csv"))
+			dfMShut = write_core_commitments(EP, inputs, FUS, :vFMSHUT, joinpath(path, "TS_FUS_maintshut.csv"))
 		end
 	end
 
 	### NON FUS CORE STARTS, SHUTS, COMMITS ###
 	if (!isempty(NONFUS) && setup["UCommit"] > 0)
-		dfNStart = write_core_behaviors(EP, inputs, :vCSTART, NONFUS, joinpath(path, "TS_NONFUS_start.csv"))
-		dfNShut = write_core_behaviors(EP, inputs, :vCSHUT, NONFUS, joinpath(path, "TS_NONFUS_shut.csv"))
-		dfNCommit = write_core_behaviors(EP, inputs, :vCCOMMIT, NONFUS, joinpath(path, "TS_NONFUS_commit.csv"))
+		dfNStart = write_core_commitments(EP, inputs, NONFUS, :vCSTART, joinpath(path, "TS_NONFUS_start.csv"))
+		dfNShut = write_core_commitments(EP, inputs, NONFUS, :vCSHUT, joinpath(path, "TS_NONFUS_shut.csv"))
+		dfNCommit = write_core_commitments(EP, inputs, NONFUS, :vCCOMMIT, joinpath(path, "TS_NONFUS_commit.csv"))
 	end
 
 	# Write dual values of certain constraints
-	write_thermal_storage_system_max_dual(path, inputs, setup, EP)
+	write_thermal_storage_system_max_dual(EP, inputs, setup, joinpath(path, "TS_System_Max_Cap_dual.csv"), msf)
+	write_thermal_storage_capacity_duals(EP, inputs, setup, joinpath(path, "TS_Capacity_Duals"), msf)
 
 end