diff --git a/.github/workflows/ClimaLandSimulations.yml b/.github/workflows/ClimaLandSimulations.yml
new file mode 100644
index 0000000000..6edcd60d35
--- /dev/null
+++ b/.github/workflows/ClimaLandSimulations.yml
@@ -0,0 +1,31 @@
+name: ClimaLandSimulations CI
+on:
+  pull_request:
+  push:
+    branches:
+      - main
+
+permissions:
+  actions: write
+  contents: read
+
+jobs:
+  lib-climalandsimulations:
+    runs-on: ubuntu-20.04
+    timeout-minutes: 45
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: '1.10'
+      - uses: julia-actions/cache@v1
+      - name: Install Julia dependencies
+        run: >
+          julia --project= -e 'using Pkg; Pkg.develop(path="$(pwd())"); Pkg.develop(path="$(pwd())/lib/ClimaLandSimulations")'
+      - name: Run the tests
+        continue-on-error: true
+        env:
+            CI_OUTPUT_DIR: output
+        run: >
+          julia --project= -e 'using Pkg; Pkg.test("ClimaLandSimulations")'
diff --git a/NEWS.md b/NEWS.md
index 62ff9235cb..af62b81546 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -23,6 +23,8 @@ v0.11.1
 -------
 - ![][badge-✨feature] Add option to profile albedo job. PR
   [#505](https://github.com/CliMA/ClimaLand.jl/pull/551)
+- ![][badge-✨feature] ClimaLandSimulations: better plots (legend and style tweaks, added cumulative ET and P), unit tests, start and end time as an optional argument. PR [#494](https://github.com/CliMA/ClimaLand.jl/pull/494)
+- ![][badge-🐛bugfix] ClimaLandSimulations: installation readme. PR [#494](https://github.com/CliMA/ClimaLand.jl/pull/494)
 
 v0.11.0
 -------
diff --git a/lib/ClimaLandSimulations/README.md b/lib/ClimaLandSimulations/README.md
index b2b2c0eef7..d24b377315 100644
--- a/lib/ClimaLandSimulations/README.md
+++ b/lib/ClimaLandSimulations/README.md
@@ -1,16 +1,27 @@
 # ClimaLandSimulations
 
 A library of methods to run ClimaLand:
-- at single sites, such as FLUXNET
+- at single sites 
 - globally (WIP)
 
 ## Installation
 
 ```jl
 julia> ] # to go into pkg mode
-pkg> add https://github.com/CliMA/ClimaLand.jl/tree/main/lib/ClimaLandSimulations # because unregistered for now
+pkg> add https://github.com/CliMA/ClimaLand.jl:lib/ClimaLandSimulations # because unregistered for now
 pkg> *backspace* # press backspace button to go back into julia mode
 julia> using ClimaLandSimulations
 ```
 
 For examples, see the [experiments](https://github.com/CliMA/ClimaLand.jl/lib/ClimaLandSimulations/experiments) folder
+
+### Development of the `ClimaLandSimulations` subpackage
+
+    cd ClimaCore/lib/ClimaLandSimulations
+
+    # Add ClimaCore to subpackage environment
+    julia --project -e 'using Pkg; Pkg.develop(path="../../")'
+
+    # Instantiate ClimaCoreMakie project environment
+    julia --project -e 'using Pkg; Pkg.instantiate()'
+    julia --project -e 'using Pkg; Pkg.test()'
diff --git a/lib/ClimaLandSimulations/experiments/ozark.jl b/lib/ClimaLandSimulations/experiments/ozark.jl
index 562e8a6f29..115ee5cca0 100644
--- a/lib/ClimaLandSimulations/experiments/ozark.jl
+++ b/lib/ClimaLandSimulations/experiments/ozark.jl
@@ -7,12 +7,20 @@ sv_test, sol_test, Y_test, p_test = run_fluxnet(
     params = ozark_default_params(; hetero_resp = hetero_resp_ozark(; b = 2)),
 )
 
+# defaults, except start time
+sv, sol, Y, p = run_fluxnet(
+    "US-MOz";
+    setup = make_setup(;
+        ozark_default_args(; t0 = Float64(125 * 3600 * 24))...,
+    ),
+)
+
 # default parameters
 sv, sol, Y, p = run_fluxnet("US-MOz")
 
 # DataFrame for input and output
 inputs, inputs_SI, inputs_commonly_used = make_inputs_df("US-MOz")
-simulation_output = make_output_df(sv, inputs)
+simulation_output = make_output_df("US-MOz", sv, inputs)
 
 # Will save figures in current repo /figures
 make_plots(inputs, simulation_output)
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_config.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_config.jl
index 0acc9a428f..3d3369c67c 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_config.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_config.jl
@@ -1,19 +1,36 @@
 export harvard_default_configs
 
 """
-    harvard_default_configs = (
-        data_link = "https://caltech.box.com/shared/static/xixaod6511cutz51ag81k1mtvy05hbol.csv",
-        local_to_UTC = 5,
-        lat = FT(42.5378), # degree
-        long = FT(-72.1715), # degree
-        atmos_h = FT(30),
+	function harvard_default_configs(;
+		data_link = "https://caltech.box.com/shared/static/xixaod6511cutz51ag81k1mtvy05hbol.csv",
+		local_to_UTC = 5,
+		lat = FT(42.5378), # degree
+		long = FT(-72.1715), # degree
+		atmos_h = FT(30),
+		)
+		return (
+			data_link = data_link,
+			local_to_UTC = local_to_UTC,
+			atmos_h = atmos_h,
+			lat = lat,
+			long = long,	
+		)
+	end
 
 Named Tuple containing default configurations for the Harvard site. 
 """
-harvard_default_configs = (
+function harvard_default_configs(;
     data_link = "https://caltech.box.com/shared/static/xixaod6511cutz51ag81k1mtvy05hbol.csv",
     local_to_UTC = 5,
     lat = FT(42.5378), # degree
     long = FT(-72.1715), # degree
     atmos_h = FT(30),
 )
+    return (
+        data_link = data_link,
+        local_to_UTC = local_to_UTC,
+        atmos_h = atmos_h,
+        lat = lat,
+        long = long,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl
index ee39585f9a..bccdbc01ef 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl
@@ -1,21 +1,33 @@
 export harvard_default_args
 
 """
-    harvard_default_args = (
-        dz_bottom = FT(1.5),
-        dz_top = FT(0.025),
-        n_stem = Int64(1),
-        n_leaf = Int64(1),
-        h_leaf = FT(12), # m
-        h_stem = FT(14), # m
-        t0 = Float64(120 * 3600 * 24),
-        dt = Float64(40),
-        n = 45,
-        )
+	function harvard_default_args(;
+		dz_bottom = FT(1.5),
+		dz_top = FT(0.025),
+		n_stem = Int64(1),
+		n_leaf = Int64(1),
+		h_leaf = FT(12), # m
+		h_stem = FT(14), # m
+		t0 = Float64(120 * 3600 * 24),
+		dt = Float64(40),
+		n = 45,
+		)
+		return (
+			dz_bottom = dz_bottom,
+			dz_top = dz_top,
+			n_stem = n_stem,
+			n_leaf = n_leaf,
+			h_stem = h_stem,
+			h_leaf = h_leaf,
+			t0 = t0,
+			dt = dt,
+			n = n,	
+		)
+	end
 
 Named Tuple containing default simulation parameter for the Harvard site. 
 """
-harvard_default_args = (
+function harvard_default_args(;
     dz_bottom = FT(1.5),
     dz_top = FT(0.025),
     n_stem = Int64(1),
@@ -26,3 +38,15 @@ harvard_default_args = (
     dt = Float64(40),
     n = 45,
 )
+    return (
+        dz_bottom = dz_bottom,
+        dz_top = dz_top,
+        n_stem = n_stem,
+        n_leaf = n_leaf,
+        h_stem = h_stem,
+        h_leaf = h_leaf,
+        t0 = t0,
+        dt = dt,
+        n = n,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_config.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_config.jl
index 7538b258f2..06ad4dd615 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_config.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_config.jl
@@ -1,20 +1,36 @@
 export ozark_default_configs
 
 """
-    ozark_default_configs = (
-        data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv",
-        local_to_UTC = 7,
-        atmos_h = FT(32),
-        lat = FT(38.7441), 
-        long = FT(-92.2000)
-        )
+	function ozark_default_configs(;
+		data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv",
+		local_to_UTC = 7,
+		atmos_h = FT(32),
+		lat = FT(38.7441),
+		long = FT(-92.2000),
+		)
+		return (
+			data_link = data_link,
+			local_to_UTC = local_to_UTC,
+			atmos_h = atmos_h,
+			lat = lat,
+			long = long,
+		)
+	end
 
 Named Tuple containing default configurations for the Ozark site. 
 """
-ozark_default_configs = (
+function ozark_default_configs(;
     data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv",
     local_to_UTC = 7,
     atmos_h = FT(32),
     lat = FT(38.7441),
     long = FT(-92.2000),
 )
+    return (
+        data_link = data_link,
+        local_to_UTC = local_to_UTC,
+        atmos_h = atmos_h,
+        lat = lat,
+        long = long,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl
index e2bb215d73..dcbff11279 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl
@@ -1,21 +1,33 @@
 export ozark_default_args
 
 """
-    ozark_default_args = (
-    dz_bottom = FT(1.5),
-    dz_top = FT(0.025),
-    n_stem = Int64(1),
-    n_leaf = Int64(1),
-    h_stem = FT(9),
-    h_leaf = FT(9.5),
-    t0 = Float64(120 * 3600 * 24), 
-    dt = Float64(120),
-    n = 15
-   )
+	function ozark_default_args(;
+		dz_bottom = FT(1.5),
+		dz_top = FT(0.025),
+		n_stem = Int64(1),
+		n_leaf = Int64(1),
+		h_stem = FT(9),
+		h_leaf = FT(9.5),
+		t0 = Float64(120 * 3600 * 24),
+		dt = Float64(120),
+		n = 15,
+		)
+		return (
+			dz_bottom = dz_bottom,
+			dz_top = dz_top,
+			n_stem = n_stem,
+			n_leaf = n_leaf,
+			h_stem = h_stem,
+			h_leaf = h_leaf,
+			t0 = t0,
+			dt = dt,
+			n = n,
+		)
+	end
 
 Named Tuple containing default simulation parameter for the Ozark site. 
 """
-ozark_default_args = (
+function ozark_default_args(;
     dz_bottom = FT(1.5),
     dz_top = FT(0.025),
     n_stem = Int64(1),
@@ -26,3 +38,15 @@ ozark_default_args = (
     dt = Float64(120),
     n = 15,
 )
+    return (
+        dz_bottom = dz_bottom,
+        dz_top = dz_top,
+        n_stem = n_stem,
+        n_leaf = n_leaf,
+        h_stem = h_stem,
+        h_leaf = h_leaf,
+        t0 = t0,
+        dt = dt,
+        n = n,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_config.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_config.jl
index 10d0989ca6..185980026f 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_config.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_config.jl
@@ -1,20 +1,36 @@
 export niwotridge_default_configs
 
 """
-    niwot_default_configs = (
-        data_link = "https://caltech.box.com/shared/static/r6gvldgabk3mvtx53gevnlaq1ztsk41i.csv",
-        local_to_UTC = 7,
-        lat = FT(40.0329), # degree,
-        long = FT(-105.5464), # degree,
-        atmos_h = FT(21.5),
-    )
+	function niwotridge_default_configs(;
+		data_link = "https://caltech.box.com/shared/static/r6gvldgabk3mvtx53gevnlaq1ztsk41i.csv",
+		local_to_UTC = 7,
+		lat = FT(40.0329), # degree,
+		long = FT(-105.5464), # degree,
+		atmos_h = FT(21.5),
+		)
+		return (
+			data_link = data_link,
+			local_to_UTC = local_to_UTC,
+			atmos_h = atmos_h,
+			lat = lat,
+			long = long,	
+		)
+	end
 
 Named Tuple containing default configurations for the Niwot Ridge site. 
 """
-niwotridge_default_configs = (
+function niwotridge_default_configs(;
     data_link = "https://caltech.box.com/shared/static/r6gvldgabk3mvtx53gevnlaq1ztsk41i.csv",
     local_to_UTC = 7,
     lat = FT(40.0329), # degree,
     long = FT(-105.5464), # degree,
     atmos_h = FT(21.5),
 )
+    return (
+        data_link = data_link,
+        local_to_UTC = local_to_UTC,
+        atmos_h = atmos_h,
+        lat = lat,
+        long = long,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl
index 3d8e853171..0d3acf86be 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl
@@ -1,21 +1,33 @@
 export niwotridge_default_args
 
 """
-    niwotridge_default_args = (
-        dz_bottom = FT(1.25),
-        dz_top = FT(0.05),
-        n_stem = Int64(1),
-        n_leaf = Int64(1),
-        h_leaf = FT(6.5), # m,
-        h_stem = FT(7.5), # m,
-        t0 = Float64(120 * 3600 * 24), # start mid year to avoid snow,
-        dt = Float64(40),
-        n = 45,
-    )
+	function niwotridge_default_args(;
+		dz_bottom = FT(1.25),
+		dz_top = FT(0.05),
+		n_stem = Int64(1),
+		n_leaf = Int64(1),
+		h_leaf = FT(6.5), # m,
+		h_stem = FT(7.5), # m,
+		t0 = Float64(120 * 3600 * 24), # start mid year to avoid snow,
+		dt = Float64(40),
+		n = 45,
+		)
+		return (
+			dz_bottom = dz_bottom,
+			dz_top = dz_top,
+			n_stem = n_stem,
+			n_leaf = n_leaf,
+			h_stem = h_stem,
+			h_leaf = h_leaf,
+			t0 = t0,
+			dt = dt,
+			n = n,	
+		)
+	end
 
 Named Tuple containing default simulation parameter for the Niwot Ridge site. 
 """
-niwotridge_default_args = (
+function niwotridge_default_args(;
     dz_bottom = FT(1.25),
     dz_top = FT(0.05),
     n_stem = Int64(1),
@@ -26,3 +38,15 @@ niwotridge_default_args = (
     dt = Float64(40),
     n = 45,
 )
+    return (
+        dz_bottom = dz_bottom,
+        dz_top = dz_top,
+        n_stem = n_stem,
+        n_leaf = n_leaf,
+        h_stem = h_stem,
+        h_leaf = h_leaf,
+        t0 = t0,
+        dt = dt,
+        n = n,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_config.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_config.jl
index d7e507b8c5..3cfd281c47 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_config.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_config.jl
@@ -1,18 +1,27 @@
 export vairaranch_default_configs
 
 """
-    vairaranch_default_configs = (
-        data_link = "https://caltech.box.com/shared/static/dx0p5idbsbrdebsda10t9pfv2lbdaz95.csv",
-        LAI_link = "https://caltech.box.com/shared/static/y5vf8s9qkoogglc1bc2eyu1k95sbjsc3.csv",
-        atmos_h = FT(2),
-        local_to_UTC = 8,
-        lat = FT(38.4133), # degree
-        long = FT(-120.9508), # degree
-    )
+function vairaranch_default_configs(;
+    data_link = "https://caltech.box.com/shared/static/dx0p5idbsbrdebsda10t9pfv2lbdaz95.csv",
+    LAI_link = "https://caltech.box.com/shared/static/y5vf8s9qkoogglc1bc2eyu1k95sbjsc3.csv",
+    atmos_h = FT(2),
+    local_to_UTC = 8,
+    lat = FT(38.4133), # degree
+    long = FT(-120.9508), # degree
+	)
+	return (
+	    data_link = data_link,
+		LAI_link = LAI_link,
+        local_to_UTC = local_to_UTC,
+        atmos_h = atmos_h,
+        lat = lat,
+        long = long,	
+	)
+end
 
 Named Tuple containing default configurations for the Vaira Ranch site. 
 """
-vairaranch_default_configs = (
+function vairaranch_default_configs(;
     data_link = "https://caltech.box.com/shared/static/dx0p5idbsbrdebsda10t9pfv2lbdaz95.csv",
     LAI_link = "https://caltech.box.com/shared/static/y5vf8s9qkoogglc1bc2eyu1k95sbjsc3.csv",
     atmos_h = FT(2),
@@ -20,3 +29,12 @@ vairaranch_default_configs = (
     lat = FT(38.4133), # degree
     long = FT(-120.9508), # degree
 )
+    return (
+        data_link = data_link,
+        LAI_link = LAI_link,
+        local_to_UTC = local_to_UTC,
+        atmos_h = atmos_h,
+        lat = lat,
+        long = long,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl
index cea1995b43..66dc5735c3 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl
@@ -1,21 +1,33 @@
 export vairaranch_default_args
 
 """
-    vairaranch_default_args = (
-        dz_bottom = FT(1.0),
-        dz_top = FT(0.05),
-        n_stem = Int64(0),
-        n_leaf = Int64(1),
-        h_leaf = FT(0.7), # m,
-        h_stem = FT(0), # m,
-        t0 = Float64(21 * 3600 * 24), # start day 21 of the year,
-        dt = Float64(40),
-        n = 45,
-    )
+	function vairaranch_default_args(;
+		dz_bottom = FT(1.0),
+		dz_top = FT(0.05),
+		n_stem = Int64(0),
+		n_leaf = Int64(1),
+		h_leaf = FT(0.7), # m,
+		h_stem = FT(0), # m,
+		t0 = Float64(21 * 3600 * 24), # start day 21 of the year,
+		dt = Float64(40),
+		n = 45,
+		)
+		return (
+			dz_bottom = dz_bottom,
+			dz_top = dz_top,
+			n_stem = n_stem,
+			n_leaf = n_leaf,
+			h_stem = h_stem,
+			h_leaf = h_leaf,
+			t0 = t0,
+			dt = dt,
+			n = n,	
+		)
+	end
 
 Named Tuple containing default simulation parameter for the Vaira Ranch site. 
 """
-vairaranch_default_args = (
+function vairaranch_default_args(;
     dz_bottom = FT(1.0),
     dz_top = FT(0.05),
     n_stem = Int64(0),
@@ -26,3 +38,15 @@ vairaranch_default_args = (
     dt = Float64(40),
     n = 45,
 )
+    return (
+        dz_bottom = dz_bottom,
+        dz_top = dz_top,
+        n_stem = n_stem,
+        n_leaf = n_leaf,
+        h_stem = h_stem,
+        h_leaf = h_leaf,
+        t0 = t0,
+        dt = dt,
+        n = n,
+    )
+end
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_config.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_config.jl
index 7330c7a794..d6b24cf3ae 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_config.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_config.jl
@@ -22,10 +22,10 @@ end
 
 function make_config(site_ID)
     default_configs = Dict(
-        "US-MOz" => ozark_default_configs,
-        "US-Ha1" => harvard_default_configs,
-        "US-NR1" => niwotridge_default_configs,
-        "US-Var" => vairaranch_default_configs,
+        "US-MOz" => ozark_default_configs(),
+        "US-Ha1" => harvard_default_configs(),
+        "US-NR1" => niwotridge_default_configs(),
+        "US-Var" => vairaranch_default_configs(),
     )
 
     if haskey(default_configs, site_ID)
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl
index 2e5bc2f1f9..5f80beee14 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl
@@ -40,10 +40,10 @@ end
 
 function make_setup(site_ID)
     default_args = Dict(
-        "US-MOz" => ozark_default_args,
-        "US-Ha1" => harvard_default_args,
-        "US-NR1" => niwotridge_default_args,
-        "US-Var" => vairaranch_default_args,
+        "US-MOz" => ozark_default_args(),
+        "US-Ha1" => harvard_default_args(),
+        "US-NR1" => niwotridge_default_args(),
+        "US-Var" => vairaranch_default_args(),
     )
 
     if haskey(default_args, site_ID)
diff --git a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl b/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl
index 460e63d319..3db7826318 100644
--- a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl
+++ b/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl
@@ -1,63 +1,85 @@
 export getoutput, make_output_df
 
 """
-    getoutput(variable::Symbol, variables::Symbol...; result = sv.saveval, depth = 1)
+    getoutput(sv, variable::Symbol, variables::Symbol...; result = sv.saveval, depth = 1)
 
 Return a vector of FT corresponding to the variable of interest at all times.
 By default, get output from sv.saveval, but user can specify e.g., result = sol.u
 By default, get surface value, but user can specify depth for e.g., soil temperature
+
+example 1: 
+julia> getoutput(sv, 1, :SW_out)
+
+example 2: 
+julia> getoutput(sv, 1, :canopy, :conductance, :gs)
+
+example 3: 
+julia> getoutput(sol, 1, :soil, :ϑ_l; result = sol.u)
 """
 function getoutput(
-    sv,
+    sv, # or sol for prognostic variables
+    depth, # 1 is surface
     variable::Symbol,
     variables::Symbol...;
-    result = sv.saveval,
-    depth = 1,
+    result = sv.saveval, # or sol.u for prognostic variables
 )
     for v in (variable, variables...)
         result = getproperty.(result, v)
     end
     return [parent(r)[depth] for r in result]
 end
-# example: getoutput(:SW_out)
-# example 2: getoutput(:canopy, :conductance, :gs)
-# example 3: getoutput(:soil, :T; result = sol.u, depth = 2)
 
 """
     make_output_df(sv, inputs)
 
 Return a dataframe containing climaland outputs
 """
-function make_output_df(sv, inputs; N_days = 60, N_spinup_days = 30)
+function make_output_df(
+    site_ID,
+    sv,
+    inputs;
+    setup = make_setup(site_ID),
+    timestepper = make_timestepper(setup),
+)
     # List of output that we want
-    output_list = (
-        (:SW_out,),
-        (:LW_out,),
-        (:canopy, :conductance, :gs),
-        (:canopy, :conductance, :transpiration),
-        (:canopy, :autotrophic_respiration, :Ra),
-        (:canopy, :photosynthesis, :GPP),
-        (:canopy, :hydraulics, :β),
-        (:canopy, :hydraulics, :area_index, :leaf),
-        (:canopy, :energy, :lhf),
-        (:soil, :turbulent_fluxes, :shf),
-        (:soil, :turbulent_fluxes, :lhf),
-        (:soil, :T),
-        (:soil, :θ_l),
-        (:soil, :turbulent_fluxes, :vapor_flux),
+    output_list = vcat(
+        (1, :SW_out),
+        (1, :LW_out),
+        (1, :canopy, :conductance, :gs),
+        (1, :canopy, :conductance, :transpiration),
+        (1, :canopy, :autotrophic_respiration, :Ra),
+        (1, :canopy, :photosynthesis, :GPP),
+        (1, :canopy, :hydraulics, :β),
+        (1, :canopy, :hydraulics, :area_index, :leaf),
+        (1, :canopy, :energy, :lhf),
+        (1, :soil, :turbulent_fluxes, :shf),
+        (1, :soil, :turbulent_fluxes, :lhf),
+        collect(map(i -> (i, :soil, :T), 1:20)), # 20 shouldn't be hard-coded, but an arg, equal to n layers
+        collect(map(i -> (i, :soil, :θ_l), 1:20)),
+        (1, :soil, :turbulent_fluxes, :vapor_flux),
     )
 
     output_vectors = [getoutput(sv, args...) for args in output_list]
     # Extract the last symbol from each tuple for column names
-    column_names = [names[end] for names in output_list]
+    column_names =
+        [Symbol(names[end], "_depth_", names[1]) for names in output_list]
+    # remove "_1"
+    for i in 1:length(column_names)
+        name = string(column_names[i])
+        if endswith(name, "_1")
+            column_names[i] = Symbol(chop(name, tail = 8))
+        end
+    end
     # Create a dictionary with simplified column names and corresponding vectors
     data_dict = Dict(zip(column_names, output_vectors))
     # Create a DataFrame from the dictionary
     climaland = DataFrame(data_dict)
     # now if I want for example GPP, I can just do df.GPP
 
-    index_t_start = 120 * 48 # we shouldn't hardcode that 120 in ozark_simulation.jl
-    index_t_end = 120 * 48 + (N_days - N_spinup_days) * 48
+    index_t_start = Int(setup.t0 / (3600 * 24) * 48)
+    index_t_end = Int(
+        index_t_start + (timestepper.N_days - timestepper.N_spinup_days) * 48,
+    )
     model_dt = inputs.DateTime[index_t_start:index_t_end]
 
     insertcols!(climaland, 1, :DateTime => model_dt)
diff --git a/lib/ClimaLandSimulations/src/utilities/make_timestepper.jl b/lib/ClimaLandSimulations/src/utilities/make_timestepper.jl
index 499a937e6e..be2316115c 100644
--- a/lib/ClimaLandSimulations/src/utilities/make_timestepper.jl
+++ b/lib/ClimaLandSimulations/src/utilities/make_timestepper.jl
@@ -30,5 +30,6 @@ function make_timestepper(
         saveat = saveat,
         timestepper = timestepper,
         ode_algo = ode_algo,
+        N_spinup_days = N_spinup_days,
     )
 end
diff --git a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl
index 82a91af4b9..5036b64d47 100644
--- a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl
+++ b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl
@@ -22,8 +22,6 @@ export timeseries_fluxes_fig,
 function timeseries_fluxes_fig(
     inputs,
     climaland,
-    index_t_start,
-    index_t_end,
     earth_param_set;
     dashboard = false,
 ) # will run for any inputs or climaland output of FLUXNET sites
@@ -32,6 +30,9 @@ function timeseries_fluxes_fig(
         WGLMakie.activate!() # for dashboards
     end
 
+    index_t_start = findfirst(isequal(climaland.DateTime[1]), inputs.DateTime)
+    index_t_end = findfirst(isequal(climaland.DateTime[end]), inputs.DateTime)
+
     fig = Figure(size = (1000, 1000)) # note: do not load Plots.jl in this branch (it is loading in plot_utils)
     fontsize_theme = Theme(fontsize = 20)
     set_theme!(fontsize_theme)
@@ -54,13 +55,13 @@ function timeseries_fluxes_fig(
         optimize_ticks(climaland.DateTime[1], climaland.DateTime[end])[1][2:(end - 1)] # first and last are weirdly placed
 
     # add plots into axis ax_C
-    p_GPP_m = scatter!(
+    p_GPP_m = lines!(
         ax_C,
         datetime2unix.(climaland.DateTime),
         climaland.GPP .* 1e6,
         color = :blue,
     )
-    p_GPP_d = scatter!(
+    p_GPP_d = lines!(
         ax_C,
         datetime2unix.(inputs.DateTime[index_t_start:index_t_end]),
         inputs.GPP[index_t_start:index_t_end] .* 1e6,
@@ -68,14 +69,14 @@ function timeseries_fluxes_fig(
     )
 
     # ax_W
-    p_ET_m = scatter!(
+    p_ET_m = lines!(
         ax_W,
         datetime2unix.(climaland.DateTime),
         (climaland.vapor_flux .* 1e3 .* 24 .* 3600) .+
         (climaland.transpiration .* 1e3 .* 24 .* 3600),
         color = :blue,
     ) # not sure about units
-    p_ET_d = scatter!(
+    p_ET_d = lines!(
         ax_W,
         datetime2unix.(inputs.DateTime[index_t_start:index_t_end]),
         inputs.LE[index_t_start:index_t_end] ./
@@ -84,13 +85,13 @@ function timeseries_fluxes_fig(
     ) # not sure units
 
     # ax_SW_OUT
-    p_SWOUT_m = scatter!(
+    p_SWOUT_m = lines!(
         ax_SWOUT,
         datetime2unix.(climaland.DateTime),
         climaland.SW_out,
         color = :blue,
     )
-    p_SWOUT_d = scatter!(
+    p_SWOUT_d = lines!(
         ax_SWOUT,
         datetime2unix.(inputs.DateTime[index_t_start:index_t_end]),
         inputs.SW_OUT[index_t_start:index_t_end],
@@ -136,8 +137,6 @@ end
 function timeseries_H2O_fig(
     inputs,
     climaland,
-    index_t_start,
-    index_t_end,
     earth_param_set;
     dashboard = false,
 ) # will run for any inputs or climaland output of FLUXNET sites
@@ -146,6 +145,10 @@ function timeseries_H2O_fig(
         WGLMakie.activate!() # for dashboards
     end
 
+    index_t_start = findfirst(isequal(climaland.DateTime[1]), inputs.DateTime)
+    index_t_end = findfirst(isequal(climaland.DateTime[end]), inputs.DateTime)
+
+
     # create an empty figure
     fig = Figure(size = (1000, 1000)) # note: do not load Plots.jl in this branch (it is loading in plot_utils)
     fontsize_theme = Theme(fontsize = 20)
@@ -175,13 +178,13 @@ function timeseries_H2O_fig(
         optimize_ticks(climaland.DateTime[1], climaland.DateTime[end])[1][2:(end - 1)] # first and last are weirdly placed
 
     # add plots into axis ax_H2O
-    p_H2O_m = scatter!(
+    p_H2O_m = lines!(
         ax_H2O,
         datetime2unix.(climaland.DateTime),
         climaland.θ_l,
         color = :green,
     )
-    p_H2O_d = scatter!(
+    p_H2O_d = lines!(
         ax_H2O,
         datetime2unix.(inputs.DateTime[index_t_start:index_t_end]),
         inputs.SWC[index_t_start:index_t_end],
@@ -197,7 +200,7 @@ function timeseries_H2O_fig(
     )
 
     # Moisture stress
-    p_MS = scatter!(
+    p_MS = lines!(
         ax_MS,
         datetime2unix.(climaland.DateTime),
         climaland.β,
@@ -205,7 +208,7 @@ function timeseries_H2O_fig(
     ) # not sure about units
 
     # Stomatal conductance
-    p_SC = scatter!(
+    p_SC = lines!(
         ax_SC,
         datetime2unix.(climaland.DateTime),
         climaland.gs,
@@ -251,19 +254,15 @@ function timeseries_H2O_fig(
 end
 
 # 3. Fingerprint plot
-function fingerprint_fig(
-    inputs,
-    climaland,
-    index_t_start,
-    index_t_end,
-    earth_param_set;
-    dashboard = false,
-) # will run for any inputs or climaland output of FLUXNET sites
+function fingerprint_fig(inputs, climaland, earth_param_set; dashboard = false) # will run for any inputs or climaland output of FLUXNET sites
 
     if dashboard == true
         WGLMakie.activate!() # for dashboards
     end
 
+    index_t_start = findfirst(isequal(climaland.DateTime[1]), inputs.DateTime)
+    index_t_end = findfirst(isequal(climaland.DateTime[end]), inputs.DateTime)
+
     fig = Figure(size = (1000, 1000))
     fontsize_theme = Theme(fontsize = 20)
     set_theme!(fontsize_theme)
@@ -275,9 +274,24 @@ function fingerprint_fig(
         xlabel = "Date",
         title = L"\text{GPP} \, (\mu\text{mol m}^{-2} \, \text{s}^{-1})",
     ) # C fluxes
-    ax_W = Axis(fig[2, 1]) # h2o fluxes
-    ax_P = Axis(fig[3, 1]) # Precip & soil moisture
-    ax_T = Axis(fig[4, 1]) # air, canopy, and soil temperature
+    ax_W = Axis(
+        fig[2, 1],
+        ylabel = "Hour of the day",
+        xlabel = "Date",
+        title = L"\text{ET (mm)}",
+    ) # h2o fluxes
+    ax_M = Axis(
+        fig[3, 1],
+        ylabel = "Hour of the day",
+        xlabel = "Date",
+        title = L"\text{soil moisture}",
+    ) # soil moisture
+    ax_R = Axis(
+        fig[4, 1],
+        ylabel = "Hour of the day",
+        xlabel = "Date",
+        title = L"\text{incoming shortwave radiation}",
+    ) # radiation
 
     # for time series, CairoMakie should allow DateTime type soon (but not yet)
     # so the 2 lines of code below are a trick to be able to use DateTime - will be removed later
@@ -293,9 +307,42 @@ function fingerprint_fig(
     )
     Colorbar(fig[1, 2], hm_GPP)
 
+    hm_ET = heatmap!(
+        ax_W,
+        datetime2unix.(DateTime.(Date.(inputs.DateTime))),
+        hour.(inputs.DateTime) .+ (minute.(inputs.DateTime) ./ 60),
+        inputs.LE ./ (LP.LH_v0(earth_param_set) * 1000) .* (1e3 * 24 * 3600),
+    )
+    Colorbar(fig[2, 2], hm_ET)
+
+    hm_M = heatmap!(
+        ax_M,
+        datetime2unix.(DateTime.(Date.(inputs.DateTime))),
+        hour.(inputs.DateTime) .+ (minute.(inputs.DateTime) ./ 60),
+        inputs.SWC,
+    )
+    Colorbar(fig[3, 2], hm_M)
+
+    hm_R = heatmap!(
+        ax_R,
+        datetime2unix.(DateTime.(Date.(inputs.DateTime))),
+        hour.(inputs.DateTime) .+ (minute.(inputs.DateTime) ./ 60),
+        inputs.SW_IN,
+    )
+    Colorbar(fig[4, 2], hm_R)
+
     ax_C.xticks[] =
         (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd"))
 
+    ax_W.xticks[] =
+        (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd"))
+
+    ax_M.xticks[] =
+        (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd"))
+
+    ax_R.xticks[] =
+        (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd"))
+
     fig
     return fig
 end
@@ -335,19 +382,16 @@ function diurnal_plot!(
     return diurnal_p
 end
 
-function diurnals_fig(
-    inputs,
-    climaland,
-    index_t_start,
-    index_t_end,
-    earth_param_set;
-    dashboard = false,
-) # will run for any inputs or climaland output of FLUXNET sites
+function diurnals_fig(inputs, climaland, earth_param_set; dashboard = false) # will run for any inputs or climaland output of FLUXNET sites
 
     if dashboard == true
         WGLMakie.activate!() # for dashboards
     end
 
+    index_t_start = findfirst(isequal(climaland.DateTime[1]), inputs.DateTime)
+    index_t_end = findfirst(isequal(climaland.DateTime[end]), inputs.DateTime)
+
+
     fig = Figure(size = (1000, 1000))
     fontsize_theme = Theme(fontsize = 20)
     set_theme!(fontsize_theme)
@@ -373,7 +417,13 @@ function diurnals_fig(
         climaland.GPP .* 1e6,
         :green,
     )
-    diurnal_plot!(fig, ax_C, climaland.DateTime, climaland.Ra .* 1e6, :black)
+    p_RA_m = diurnal_plot!(
+        fig,
+        ax_C,
+        climaland.DateTime,
+        climaland.Ra .* 1e6,
+        :black,
+    )
     # data
     p_GPP_d = diurnal_plot!(
         fig,
@@ -426,27 +476,27 @@ function diurnals_fig(
 
     axislegend(
         ax_C,
-        [p_GPP_d, p_GPP_m],
-        ["Observations", "ClimaLand"],
+        [p_GPP_d, p_GPP_m, p_RA_m],
+        ["GPP Obs.", "GPP model", "Ra model"],
         "",
         position = :rt,
-        orientation = :horizontal,
+        orientation = :vertical,
     )
     axislegend(
         ax_W,
         [p_ET_d, p_ET_m],
-        ["Observations", "ClimaLand"],
+        ["ET obs.", "ET model"],
         "",
         position = :rt,
-        orientation = :horizontal,
+        orientation = :vertical,
     )
     axislegend(
         ax_E,
         [p_SWout_d, p_SWout_m],
-        ["Observations", "ClimaLand"],
+        ["SWout obs.", "SWout model"],
         "",
         position = :rt,
-        orientation = :horizontal,
+        orientation = :vertical,
     )
 
     hidexdecorations!(ax_C)
@@ -457,6 +507,49 @@ function diurnals_fig(
 end
 
 # 5. Cumulative P and ET
+function cumulative_H2O_fig(
+    inputs,
+    climaland,
+    earth_param_set;
+    dashboard = false,
+)
+
+    if dashboard == true
+        WGLMakie.activate!() # for dashboards
+    end
+
+    index_t_start = findfirst(isequal(climaland.DateTime[1]), inputs.DateTime)
+    index_t_end = findfirst(isequal(climaland.DateTime[end]), inputs.DateTime)
+
+
+    fig = Figure(size = (1000, 1000))
+    fontsize_theme = Theme(fontsize = 20)
+    set_theme!(fontsize_theme)
+
+    ax = Axis(fig[1, 1], ylabel = "Cumulative water flux (mm)")
+
+    ET_m = climaland.transpiration .* 1e3 .* 24 .* 3600
+    ET_obs =
+        inputs.LE[index_t_start:index_t_end] ./
+        (LP.LH_v0(earth_param_set) * 1000) .* (1e3 * 24 * 3600)
+    P_obs = inputs.P[index_t_start:index_t_end] .* 1e3 .* 24 .* 3600
+
+    p_P = lines!(ax, (1 / 48):(1 / 48):(length(P_obs) / 48), cumsum(P_obs)) # Precip
+    p_ET_m = lines!(ax, (1 / 48):(1 / 48):(length(ET_m) / 48), cumsum(ET_m)) # ET model
+    p_ET_d = lines!(ax, (1 / 48):(1 / 48):(length(ET_obs) / 48), cumsum(ET_obs)) # ET observations
+
+    axislegend(
+        ax,
+        [p_P, p_ET_m, p_ET_d],
+        ["Precipitation", "ET modeled", "ET observed"],
+        "",
+        position = :lt,
+        orientation = :vertical,
+    )
+
+    fig
+    return fig
+end
 
 # 6. Energy balance closure (L + H = Rn - G)
 
@@ -475,38 +568,13 @@ function make_plots(
     end
 
     earth_param_set = LP.LandParameters(FT)
-    # below shouldn't be hardcoded!
-    index_t_start = 120 * 48 # we shouldn't hardcode that 120 in ozark_simulation.jl
-    index_t_end = 120 * 48 + (60 - 30) * 48
-
-    fig1 = timeseries_fluxes_fig(
-        inputs,
-        climaland,
-        index_t_start,
-        index_t_end,
-        earth_param_set,
-    )
-    fig2 = timeseries_H2O_fig(
-        inputs,
-        climaland,
-        index_t_start,
-        index_t_end,
-        earth_param_set,
-    )
-    fig3 = fingerprint_fig(
-        inputs,
-        climaland,
-        index_t_start,
-        index_t_end,
-        earth_param_set,
-    )
-    fig4 = diurnals_fig(
-        inputs,
-        climaland,
-        index_t_start,
-        index_t_end,
-        earth_param_set,
-    )
+    args = (inputs, climaland, earth_param_set)
+
+    fig1 = timeseries_fluxes_fig(args...)
+    fig2 = timeseries_H2O_fig(args...)
+    fig3 = fingerprint_fig(args...)
+    fig4 = diurnals_fig(args...)
+    fig5 = cumulative_H2O_fig(args...)
 
     if save_fig == true
         if isdir("figures")
@@ -520,11 +588,12 @@ function make_plots(
             "timeseries_H2O.pdf",
             "fingerprint.pdf",
             "diurnals.pdf",
+            "cumulative_water.pdf",
         ]
 
         [
             save(joinpath("figures", name), fig) for
-            (name, fig) in zip(names, [fig1, fig2, fig3, fig4])
+            (name, fig) in zip(names, [fig1, fig2, fig3, fig4, fig5])
         ]
 
         return nothing
@@ -536,6 +605,7 @@ function make_plots(
             water = fig2,
             fingerprint = fig3,
             diurnals = fig4,
+            cumulative_water = fig5,
         )
     end
 
diff --git a/lib/ClimaLandSimulations/test/Project.toml b/lib/ClimaLandSimulations/test/Project.toml
new file mode 100644
index 0000000000..1d19e585e2
--- /dev/null
+++ b/lib/ClimaLandSimulations/test/Project.toml
@@ -0,0 +1,3 @@
+[deps]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ClimaLandSimulations = "348a0bd3-1299-4261-8002-d2cd97df6055"
diff --git a/lib/ClimaLandSimulations/test/runtests.jl b/lib/ClimaLandSimulations/test/runtests.jl
index 6184309d28..655115dee5 100644
--- a/lib/ClimaLandSimulations/test/runtests.jl
+++ b/lib/ClimaLandSimulations/test/runtests.jl
@@ -1,3 +1,31 @@
-using ClimaLand # rename to ClimaLand WIP
-using ClimaLandSite
 using Test
+using ClimaLandSimulations
+using ClimaLandSimulations.Fluxnet
+
+@testset "Fluxnet single site" begin
+    sv, sol, Y, p = run_fluxnet("US-MOz")
+    @test typeof(sv) ==
+          @NamedTuple{t::Vector{Float64}, saveval::Vector{NamedTuple}}
+end
+
+@testset "Fluxnet single site, custom param" begin
+    sv, sol, Y, p = run_fluxnet(
+        "US-MOz";
+        params = ozark_default_params(;
+            hetero_resp = hetero_resp_ozark(; b = 2),
+        ),
+    )
+    @test typeof(sv) ==
+          @NamedTuple{t::Vector{Float64}, saveval::Vector{NamedTuple}}
+end
+
+@testset "Fluxnet single site, custom start time" begin
+    sv, sol, Y, p = run_fluxnet(
+        "US-MOz";
+        setup = make_setup(;
+            ozark_default_args(; t0 = Float64(125 * 3600 * 24))...,
+        ),
+    )
+    @test typeof(sv) ==
+          @NamedTuple{t::Vector{Float64}, saveval::Vector{NamedTuple}}
+end