diff --git a/scripts/runcode.jl b/scripts/runcode.jl
index d4fdeb9..6821e1f 100644
--- a/scripts/runcode.jl
+++ b/scripts/runcode.jl
@@ -224,7 +224,6 @@ julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialCo
             # initialize vector to store quantities
             param_saved_var.storeTime = Vector{Float64}([0])
             param_saved_var.storeTemp = Vector{Float64}([T_0])
-            param_saved_var.storeEps_x = Vector{Float64}([eps_x0])
 
             @info("IC_Finder parameters: $(param_IC_Finder)")
 
@@ -266,14 +265,13 @@ julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialCo
         print_timer_log(io, to)
         close(io)
         df = DataFrame(sol)
-        write_csv(df, erupt_saved, path)
+        write_csv(df, erupt_saved, path, composition)
         if plotfig
             title = latexstringtitle(composition, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth)
             plot_combined_fig(df, path, title)
-            plot_dual_axis(df, param_saved_var.storeTime, param_saved_var.storeEps_x, path)
+            # plot_figs(df, path)
+            plot_dual_axis(df, path)
         end
-        # write_ϵx_csv(param_saved_var.storeTime, param_saved_var.storeEps_x, path)
-        # plot_ϵx(param_saved_var.storeTime, param_saved_var.storeEps_x, path)
         return df
     end
     return df
diff --git a/src/Chamber.jl b/src/Chamber.jl
index b0be27e..61c02d5 100644
--- a/src/Chamber.jl
+++ b/src/Chamber.jl
@@ -140,8 +140,6 @@ export @memoize,
     latexstringtitle,
     plot_combined_fig,
     plot_dual_axis,
-    plot_ϵx,
-    write_ϵx_csv,
     check_for_duplicates,
     chamber
 
diff --git a/src/plot_figs.jl b/src/plot_figs.jl
index d4ded76..212d771 100644
--- a/src/plot_figs.jl
+++ b/src/plot_figs.jl
@@ -5,6 +5,7 @@ function plot_sol(df, x_axis_col, y_axis_col)
         xaxis=x_axis_col,
         yaxis=y_axis_col,
         linewidth=2,
+        legend=false,
         marker=(:x, 3, Plots.stroke(2)),
     )
 end
@@ -21,6 +22,7 @@ function plot_sol_year_unit(df, xlabel, ylabel)
         "time" => "Time (yr)",
         "P+dP" => "Pressure (MPa)",
         "T" => "Temperature (K)",
+        "eps_g" => L"𝜺_g",
         "V" => "Volume (km³)",
         "total_mass" => "Total mass (×10¹² kg)",
     )
@@ -47,27 +49,50 @@ end
 
 function plot_combined_fig(df::DataFrame, path::String, title::LaTeXString)::Nothing
     ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
+    # plot_names = ["P+dP", "T", "eps_g", "V", "X_CO2", "total_mass"]
+    # plts = [plot_sol_year_unit(df, "time", l) for l in plot_names]
+    # combined_plot = plot(plts..., layout=(2, 3), size=(900, 500), dpi=200, margin = 3Plots.mm, plot_title=title, plot_titlefontsize=13)
+    # # add eps_x in eps_g plot
+    # color2 = RGB(255/255, 100/255, 0/255)
+    # twin_subplot3 = twinx(combined_plot[3])
+    # plot!(twin_subplot3, df[:, 1]/31536000, df[!, 12];linewidth=2, legend=false, linecolor=color2,fg_color_guide=color2,fg_color_text=color2, yaxis=L"𝜺_X")
+    # savefig(combined_plot, joinpath(path, "combined_plots.png"))
+
+    ###
+    combined_plot = plot(layout=(2, 3), size=(900, 500), dpi=400, margin=3Plots.mm, plot_title=title, plot_titlefontsize=13, legend=false)
     plot_names = ["P+dP", "T", "eps_g", "V", "X_CO2", "total_mass"]
-    plts = [plot_sol_year_unit(df, "time", l) for l in plot_names]
-    combined_plot = plot(plts..., layout=(2, 3), size=(900, 500), dpi=200, margin = 3Plots.mm, plot_title=title, plot_titlefontsize=13)
+    # Units of outputs: Time -> sec, Presure -> Pa, Volume -> m³
+    divisors = Dict(
+        "time" => 31536000,
+        "P+dP" => 1e6,
+        "V" => 1e9,
+        "total_mass" => 1e12,
+    )
+    labels = Dict(
+        "time" => "Time (yr)",
+        "P+dP" => "Pressure (MPa)",
+        "T" => "Temperature (K)",
+        "eps_g" => L"𝜺_g",
+        "V" => "Volume (km³)",
+        "total_mass" => "Total mass (×10¹² kg)",
+    )
+    convert_unit(label) = df[:, label] / get(divisors, label, 1.0)
+    getaxis(label) = get(labels, label, label)
+    x_data = convert_unit("time")
+    for (i, label) in enumerate(plot_names)
+        subplot = combined_plot[i]
+        y_data = convert_unit(label)
+        plot!(subplot, x_data, y_data; xaxis=getaxis("time"), yaxis=getaxis(label), linewidth=2, legend=false)
+    end
+    # add eps_x in eps_g plot
+    color2 = RGB(255/255, 100/255, 0/255)
+    twin_subplot3 = twinx(combined_plot[3])
+    plot!(twin_subplot3, x_data, df[!, "eps_x"]; linewidth=2, yaxis=L"𝜺_X", legend=false, linecolor=color2, fg_color_guide=color2, fg_color_text=color2)
     savefig(combined_plot, joinpath(path, "combined_plots.png"))
     return nothing
 end
 
-function plot_ϵx(storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing
-    ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
-    savefig(plot(
-            storeTime, storeEps_x;
-            xaxis="time",
-            yaxis="eps_x",
-            label="Chamber.jl",
-            linewidth=2,
-            marker=(:x, 3, Plots.stroke(2))),
-        "$path/eps_x.png")
-    return nothing
-end
-
-function plot_dual_axis(df::DataFrame, storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing
+function plot_dual_axis(df::DataFrame, path::String)::Nothing
     color1 = RGB(20/255, 125/255, 204/255)
     color2 = RGB(255/255, 100/255, 0/255)
     ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
@@ -115,8 +140,8 @@ function plot_dual_axis(df::DataFrame, storeTime::Vector{Float64}, storeEps_x::V
         )
     p2 = plot!(
         twinx(),
-        storeTime/31536000,
-        storeEps_x;
+        df[!, 1]/31536000,
+        df[!, 12];
         common_opt2...,
         yaxis=L"𝜺_X",
     )
diff --git a/src/runcode-func.jl b/src/runcode-func.jl
index a706440..cd41dc8 100644
--- a/src/runcode-func.jl
+++ b/src/runcode-func.jl
@@ -31,7 +31,6 @@ function odeChamber(
     composition = param.composition
     storeTime = param_saved_var.storeTime
     storeTemp = param_saved_var.storeTemp
-    storeEps_x = param_saved_var.storeEps_x
     phase = param_saved_var.phase
     c_x, c_m = param.c_x, param.c_m
     L_e, L_m = param.L_e, param.L_m
@@ -67,18 +66,15 @@ function odeChamber(
     =#
     if storeTime[end] == t
         storeTemp[end] = T
-        storeEps_x[end] = eps_x
     elseif t != 0
         storeTime = [storeTime; t]
         storeTemp = [storeTemp; T]
-        storeEps_x = [storeEps_x; eps_x]
     end
 
     cross = findfirst(!=(0), diff(sign.(diff(storeTime))))
     if cross !== nothing
         cross_time = storeTime[end]
         storeTemp = [storeTemp[storeTime .< cross_time]; storeTemp[end]]
-        storeEps_x = [storeEps_x[storeTime .< cross_time]; storeEps_x[end]]
         storeTime = [storeTime[storeTime .< cross_time]; cross_time]
     end
 
@@ -204,7 +200,6 @@ function odeChamber(
     du[10] = Mdot_c_in - Mdot_c_out
     param_saved_var.storeTime = storeTime
     param_saved_var.storeTemp = storeTemp
-    param_saved_var.storeEps_x = storeEps_x
     return du
 end
 
@@ -312,7 +307,6 @@ function affect!(
 )
     composition = param.composition
     param_saved_var.storeTemp = param_saved_var.storeTemp[param_saved_var.storeTime .< int.t]
-    param_saved_var.storeEps_x = param_saved_var.storeEps_x[param_saved_var.storeTime .< int.t]
     param_saved_var.storeTime = param_saved_var.storeTime[param_saved_var.storeTime .< int.t]
 
     if param.dP_lit_dt_0 == 0
@@ -342,6 +336,7 @@ function affect!(
         record_erupt_end(int.t, erupt_saved, param)
         @info("*event idx: $idx\n  time: $(int.t), Finished an eruption...")
     elseif idx == 6 || idx == 8
+        println("idx = 8..")
         phase_here = param_saved_var.phase
         @info(
             "*event idx: $idx\n  time: $(int.t), starting ic finder for conversion of phase, phase_here: $phase_here",
diff --git a/src/utils.jl b/src/utils.jl
index d82d1ad..f0a52bb 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -137,7 +137,6 @@ end
     storeSumk_2::Vector{T} = []
     storeSumk_old::Vector{T} = []
     storeSumk_2_old::Vector{T} = []
-    storeEps_x::Vector{T} = []
 end
 
 @with_kw mutable struct EruptSaved{T}
diff --git a/src/write_csv.jl b/src/write_csv.jl
index 87c56a7..ae3197b 100644
--- a/src/write_csv.jl
+++ b/src/write_csv.jl
@@ -1,4 +1,4 @@
-function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String)::Nothing
+function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String, composition::Union{Silicic,Mafic})::Nothing
     number_of_data = size(df, 1)
     println("number_of_data: $number_of_data")
     names = [
@@ -15,6 +15,12 @@ function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String
     ]
     rename!(df, ["timestamp" => "time"])
     rename!(df, ["value$i" => names[i] for i in 1:10])
+    # add eps_x column
+    m_h2o = df[:, "total_mass_H2O"] ./ df[:, "total_mass"]
+    m_co2 = df[:, "total_mass_CO2"] ./ df[:, "total_mass"]
+    crystal_fraction_eps_x′(T, P, mH2O, mCO2) = crystal_fraction_eps_x(composition, T, P, mH2O, mCO2)
+    df[!, "eps_x"] = crystal_fraction_eps_x′.(df[:, "T"], df[:, "P+dP"], m_h2o, m_co2)
+
     n = length(erupt_saved.time)
     df_erupt = DataFrame(
         time_of_eruption = erupt_saved.time,
@@ -26,9 +32,3 @@ function write_csv(df::DataFrame, erupt_saved::EruptSaved{Float64}, path::String
     CSV.write("$path/eruptions.csv", df_erupt)
     return nothing
 end
-
-function write_ϵx_csv(storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing
-    df = DataFrame(Time = storeTime, eps_x = storeEps_x)
-    CSV.write("$path/eps_x.csv", df)
-    return nothing
-end