diff --git a/.buildkite/longruns_gpu/pipeline.yml b/.buildkite/longruns_gpu/pipeline.yml
index f023ff297f..0082fbac2f 100644
--- a/.buildkite/longruns_gpu/pipeline.yml
+++ b/.buildkite/longruns_gpu/pipeline.yml
@@ -39,6 +39,7 @@ steps:
 
   - group: "Global Land Models"
     if: build.env("LONGER_RUN") == null
+    if: build.env("PERTURBED_RUN") == null
     steps:
 
       - label: ":snow_capped_mountain: Snowy Land"
@@ -52,17 +53,6 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
-      - label: ":snow_capped_mountain: Low-Res Snowy Land"
-        command:
-          - julia --color=yes --project=.buildkite experiments/long_runs/low_res_snowy_land.jl
-        artifact_paths:
-          - "lowres_snowy_land_longrun_gpu/*png"
-        agents:
-          slurm_gpus: 1
-          slurm_time: 3:00:00
-        env:
-          CLIMACOMMS_DEVICE: "CUDA"
-
       - label: ":sunglasses: California regional simulation"
         command:
           - julia --color=yes --project=.buildkite experiments/long_runs/land_region.jl
@@ -93,6 +83,34 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
+  - group: "Perturbed Low-res Global Land Models"
+    if: build.env("PERTURBED_RUN") != null
+    steps:
+
+      - label: ":snow_capped_mountain: Low-Res Snowy Land Perturbed Temperature"
+        command:
+          - julia --color=yes --project=.buildkite experiments/long_runs/low_res_snowy_land_temp.jl
+        artifact_paths:
+          - "lowres_snowy_land_longrun_temp_gpu/*png"
+        agents:
+          slurm_gpus: 1
+          slurm_time: 3:00:00
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+          PERTURBED_RUN: ""
+ 
+      - label: ":snow_capped_mountain: Low-Res Snowy Land Perturbed RH"
+        command:
+          - julia --color=yes --project=.buildkite experiments/long_runs/low_res_snowy_land_rh.jl
+        artifact_paths:
+          - "lowres_snowy_land_longrun_rh_gpu/*png"
+        agents:
+          slurm_gpus: 1
+          slurm_time: 3:00:00
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+          PERTURBED_RUN: ""
+
   - group: "Longer runs of Global Land Models"
     if: build.env("LONGER_RUN") != null
     steps:
diff --git a/experiments/long_runs/low_res_snowy_land.jl b/experiments/long_runs/low_res_snowy_land_rh.jl
similarity index 97%
rename from experiments/long_runs/low_res_snowy_land.jl
rename to experiments/long_runs/low_res_snowy_land_rh.jl
index e093ef8c56..761d22e27c 100644
--- a/experiments/long_runs/low_res_snowy_land.jl
+++ b/experiments/long_runs/low_res_snowy_land_rh.jl
@@ -45,7 +45,7 @@ context = ClimaComms.context()
 ClimaComms.init(context)
 device = ClimaComms.device()
 device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
-root_path = "lowres_snowy_land_longrun_$(device_suffix)"
+root_path = "lowres_snowy_land_longrun_rh_$(device_suffix)"
 diagnostics_outdir = joinpath(root_path, "global_diagnostics")
 outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
@@ -58,12 +58,12 @@ function setup_model(FT, start_date, stop_date, Δt, domain, earth_param_set)
         context,
         lowres = true,
     )
-    atmos, radiation = ClimaLand.prescribed_perturbed_forcing_era5(
+    atmos, radiation = ClimaLand.prescribed_perturbed_rh_era5(
         era5_ncdata_path,
         surface_space,
         start_date,
         earth_param_set,
-        5.0,
+        -0.2,
         FT;
         max_wind_speed = 25.0,
         time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
diff --git a/experiments/long_runs/low_res_snowy_land_temp.jl b/experiments/long_runs/low_res_snowy_land_temp.jl
new file mode 100644
index 0000000000..3047b421ad
--- /dev/null
+++ b/experiments/long_runs/low_res_snowy_land_temp.jl
@@ -0,0 +1,185 @@
+# # Global run of land model at low resolution
+
+# The code sets up and runs ClimaLand v1, which
+# includes soil, canopy, and snow, on a spherical domain,
+# using low resolution ERA5 data as forcing.
+
+# Simulation Setup
+# Number of spatial elements: 30 in horizontal, 15 in vertical
+# Soil depth: 50 m
+# Simulation duration: 730 d
+# Timestep: 450 s
+# Timestepper: ARS111
+# Fixed number of iterations: 3
+# Jacobian update: every new Newton iteration
+# Atmos forcing update: every 3 hours
+
+import ClimaComms
+ClimaComms.@import_required_backends
+using ClimaUtilities.ClimaArtifacts
+import ClimaUtilities.TimeManager: ITime, date
+
+import ClimaDiagnostics
+import ClimaUtilities
+
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
+import ClimaUtilities.ClimaArtifacts: @clima_artifact
+import ClimaParams as CP
+using ClimaCore
+using ClimaLand
+using ClimaLand.Snow
+using ClimaLand.Soil
+using ClimaLand.Canopy
+import ClimaLand
+import ClimaLand.Parameters as LP
+import ClimaLand.Simulations: LandSimulation, solve!
+
+using Dates
+
+using CairoMakie, GeoMakie, Poppler_jll, ClimaAnalysis
+import ClimaLand.LandSimVis as LandSimVis
+
+const FT = Float64;
+context = ClimaComms.context()
+ClimaComms.init(context)
+device = ClimaComms.device()
+device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
+root_path = "lowres_snowy_land_longrun_temp_$(device_suffix)"
+diagnostics_outdir = joinpath(root_path, "global_diagnostics")
+outdir =
+    ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
+
+function setup_model(FT, start_date, stop_date, Δt, domain, earth_param_set)
+    surface_domain = ClimaLand.Domains.obtain_surface_domain(domain)
+    surface_space = domain.space.surface
+    # Forcing data
+    era5_ncdata_path = ClimaLand.Artifacts.era5_land_forcing_data2008_path(;
+        context,
+        lowres = true,
+    )
+    atmos, radiation = ClimaLand.prescribed_perturbed_temperature_era5(
+        era5_ncdata_path,
+        surface_space,
+        start_date,
+        earth_param_set,
+        10.0,
+        FT;
+        max_wind_speed = 25.0,
+        time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+    )
+    forcing = (; atmos, radiation)
+
+    # Read in LAI from MODIS data
+    modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(;
+        context = nothing,
+        start_date,
+        end_date = stop_date,
+    )
+    LAI = ClimaLand.prescribed_lai_modis(
+        modis_lai_ncdata_path,
+        surface_space,
+        start_date;
+        time_interpolation_method = LinearInterpolation(),
+    )
+
+    # Overwrite some defaults for the canopy model
+    # Energy model
+    ac_canopy = FT(2.5e3)
+    energy = Canopy.BigLeafEnergyModel{FT}(; ac_canopy)
+
+    # Plant hydraulics
+    a = FT(0.2 * 0.0098) # 1/m
+    retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
+    hydraulics =
+        Canopy.PlantHydraulicsModel{FT}(surface_domain, LAI; retention_model)
+
+    # Roughness lengths
+    h_canopy = hydraulics.compartment_surfaces[end]
+    z0_m = FT(0.13) * h_canopy
+    z0_b = FT(0.1) * z0_m
+
+    ground = ClimaLand.PrognosticGroundConditions{FT}()
+    canopy_forcing = (; atmos, radiation, ground)
+
+    canopy = ClimaLand.Canopy.CanopyModel{FT}(
+        surface_domain,
+        canopy_forcing,
+        LAI,
+        earth_param_set;
+        prognostic_land_components = (:canopy, :snow, :soil, :soilco2),
+        energy,
+        hydraulics,
+        z_0m = z0_m,
+        z_0b = z0_b,
+    )
+
+    # Snow model setup
+    # Set β = 0 in order to regain model without density dependence
+    α_snow = Snow.ZenithAngleAlbedoModel(
+        FT(0.64),
+        FT(0.06),
+        FT(2);
+        β = FT(0.4),
+        x0 = FT(0.2),
+    )
+    horz_degree_res =
+        sum(ClimaLand.Domains.average_horizontal_resolution_degrees(domain)) / 2 # mean of resolution in latitude and longitude, in degrees
+    scf = Snow.WuWuSnowCoverFractionModel(
+        FT(0.08),
+        FT(1.77),
+        FT(1.0),
+        horz_degree_res,
+    )
+    snow = Snow.SnowModel(
+        FT,
+        surface_domain,
+        forcing,
+        earth_param_set,
+        Δt;
+        prognostic_land_components = (:canopy, :snow, :soil, :soilco2),
+        α_snow,
+        scf,
+    )
+
+    # Construct the land model with all default components except for snow
+    land =
+        LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; snow, canopy)
+    return land
+end
+# Note that since the Northern hemisphere's winter season is defined as DJF,
+# we simulate from and until the beginning of
+# March so that a full season is included in seasonal metrics.
+start_date = DateTime("2008-03-01")
+stop_date = DateTime("2010-03-02")
+Δt = 450.0
+nelements = (30, 15)
+domain = ClimaLand.Domains.global_domain(
+    FT;
+    context,
+    nelements,
+    mask_threshold = FT(0.99),
+)
+params = LP.LandParameters(FT)
+model = setup_model(FT, start_date, stop_date, Δt, domain, params)
+user_callbacks = (
+    ClimaLand.NaNCheckCallback(
+        Dates.Month(6),
+        start_date,
+        ITime(Δt, epoch = start_date),
+        mask = ClimaLand.Domains.landsea_mask(ClimaLand.get_domain(model)),
+    ),
+    ClimaLand.ReportCallback(10000),
+)
+simulation =
+    LandSimulation(start_date, stop_date, Δt, model; user_callbacks, outdir)
+@info "Run: Global Soil-Canopy-Snow Model"
+@info "Resolution: $nelements"
+@info "Timestep: $Δt s"
+@info "Start Date: $start_date"
+@info "Stop Date: $stop_date"
+ClimaLand.Simulations.solve!(simulation)
+
+LandSimVis.make_annual_timeseries(simulation; savedir = root_path)
+LandSimVis.make_heatmaps(simulation; savedir = root_path, date = stop_date)
+LandSimVis.make_leaderboard_plots(simulation; savedir = root_path)
diff --git a/src/ClimaLand.jl b/src/ClimaLand.jl
index 7aa8678a6d..281a180229 100644
--- a/src/ClimaLand.jl
+++ b/src/ClimaLand.jl
@@ -24,6 +24,7 @@ include("shared_utilities/checkpoints.jl")
 include("shared_utilities/utils.jl")
 include("shared_utilities/models.jl")
 include("shared_utilities/drivers.jl")
+include("shared_utilities/perturbed_drivers.jl")
 include("shared_utilities/boundary_conditions.jl")
 include("shared_utilities/sources.jl")
 include("shared_utilities/implicit_timestepping.jl")
diff --git a/src/shared_utilities/drivers.jl b/src/shared_utilities/drivers.jl
index 89a762f65d..6d2f42f845 100644
--- a/src/shared_utilities/drivers.jl
+++ b/src/shared_utilities/drivers.jl
@@ -37,7 +37,8 @@ export AbstractAtmosphericDrivers,
     make_update_drivers,
     prescribed_lai_era5,
     prescribed_forcing_era5,
-    prescribed_perturbed_forcing_era5,
+    prescribed_perturbed_temperature_era5,
+    prescribed_perturbed_rh_era5,
     prescribed_analytic_forcing,
     default_zenith_angle
 
@@ -1045,47 +1046,6 @@ function specific_humidity_from_dewpoint(
     return q
 end
 
-
-"""
-    perturbed_specific_humidity_from_dewpoint(dewpoint_temperature, temperature, air_pressure, earth_param_set, ΔT)
-
-Estimates the perturbed specific humidity given the true dewpoint temperature, true temperature of the air
-in Kelvin, and true air pressure in Pa, along with a perturbation to the temperature ΔT, and the ClimaLand
-`earth_param_set`. We first compute relative humidity, and then using RH and the perturbed air temperature,
-equal to  true temperature+ΔT, compute a perturbed `q`.
-
-Relative humidity is computed using the Magnus formula.
-
-For more information on the Magnus formula, see e.g.
-Lawrence, Mark G. (1 February 2005).
-"The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air:
-A Simple Conversion and Applications".
-Bulletin of the American Meteorological Society. 86 (2): 225–234.
-"""
-function perturbed_specific_humidity_from_dewpoint(
-    T_dew_air::data_FT,
-    T_air::data_FT,
-    P_air::data_FT,
-    earth_param_set,
-    ΔT,
-) where {data_FT <: Real}
-    thermo_params = LP.thermodynamic_parameters(earth_param_set)
-    _T_freeze = LP.T_freeze(earth_param_set)
-    sim_FT = typeof(_T_freeze)
-    # Obtain the relative humidity. This function requires temperatures in Celsius
-    rh::sim_FT = rh_from_dewpoint(
-        sim_FT(T_dew_air) - _T_freeze,
-        sim_FT(T_air) - _T_freeze,
-    )
-    q = Thermodynamics.q_vap_from_RH_liquid(
-        thermo_params,
-        sim_FT(P_air),
-        sim_FT(T_air + ΔT),
-        rh,
-    )
-    return q
-end
-
 """
     rh_from_dewpoint(Td_C, T_C)
 
@@ -1687,195 +1647,6 @@ function prescribed_forcing_era5(
     return (; atmos, radiation)
 end
 
-"""
-     prescribed_perturbed_forcing_era5(era5_ncdata_path,
-                             surface_space,
-                             start_date,
-                             earth_param_set,
-                             ΔT,
-                             FT;
-                             gustiness=1,
-                             max_wind_speed = nothing,
-                             c_co2 = TimeVaryingInput((t) -> 4.2e-4),
-                             time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
-                             regridder_type = :InterpolationsRegridder,
-                             interpolation_method = Interpolations.Constant(),)
-
-A helper function which constructs the `PrescribedAtmosphere` and `PrescribedRadiativeFluxes`
-from a file path pointing to the ERA5 data in a netcdf file, the surface_space, the start date,
-and the earth_param_set, applying a change in the instantaneous temperature at each point
-of ΔT, while keeping the relative humidity fixed.
-
-Please see the documentation for `prescribed_forcing_era5` for more information.
-"""
-function prescribed_perturbed_forcing_era5(
-    era5_ncdata_path,
-    surface_space,
-    start_date,
-    earth_param_set,
-    ΔT,
-    FT;
-    gustiness = 1,
-    max_wind_speed = nothing,
-    c_co2 = TimeVaryingInput((t) -> 4.2e-4),
-    time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
-    regridder_type = :InterpolationsRegridder,
-    interpolation_method = Interpolations.Constant(),
-)
-    # Pass a list of files in all cases
-    era5_ncdata_path isa String && (era5_ncdata_path = [era5_ncdata_path])
-    _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
-    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
-    precip = TimeVaryingInput(
-        [era5_ncdata_path, era5_ncdata_path],
-        ["mtpr", "msr"],
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
-        method = time_interpolation_method,
-        compose_function = (mtpr, msr) -> min.(mtpr .- msr, Float32(0)),
-    )
-    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
-    snow_precip = TimeVaryingInput(
-        era5_ncdata_path,
-        "msr",
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
-        method = time_interpolation_method,
-    )
-    if max_wind_speed isa Nothing
-        wind_function = (u, v) -> sqrt.(u .^ 2 .+ v .^ 2)
-    else
-        wind_function = (u, v) -> min.(sqrt.(u .^ 2 .+ v .^ 2), max_wind_speed)
-    end
-
-    u_atmos = TimeVaryingInput(
-        [era5_ncdata_path, era5_ncdata_path],
-        ["u10", "v10"],
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        compose_function = wind_function,
-        method = time_interpolation_method,
-    )
-    specific_humidity(Td, T, P; params = earth_param_set, ΔT = ΔT) =
-        ClimaLand.perturbed_specific_humidity_from_dewpoint.(
-            Td,
-            T,
-            P,
-            params,
-            ΔT,
-        )
-    q_atmos = TimeVaryingInput(
-        [era5_ncdata_path, era5_ncdata_path, era5_ncdata_path],
-        ["d2m", "t2m", "sp"],
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        compose_function = specific_humidity,
-        method = time_interpolation_method,
-    )
-    P_atmos = TimeVaryingInput(
-        era5_ncdata_path,
-        "sp",
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        method = time_interpolation_method,
-    )
-    perturb_temp(T; ΔT = ΔT) = T .+ ΔT
-    T_atmos = TimeVaryingInput(
-        era5_ncdata_path,
-        "t2m",
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        file_reader_kwargs = (; preprocess_func = perturb_temp),
-        method = time_interpolation_method,
-    )
-    h_atmos = FT(10)
-
-    atmos = PrescribedAtmosphere(
-        precip,
-        snow_precip,
-        T_atmos,
-        u_atmos,
-        q_atmos,
-        P_atmos,
-        start_date,
-        h_atmos,
-        earth_param_set;
-        gustiness = FT(gustiness),
-        c_co2 = c_co2,
-    )
-
-    SW_d = TimeVaryingInput(
-        era5_ncdata_path,
-        "msdwswrf",
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        method = time_interpolation_method,
-    )
-    function compute_diffuse_fraction(total, direct)
-        diff = max(total - direct, Float32(0))
-        return min(diff / (total + eps(Float32)), Float32(1))
-    end
-    function compute_diffuse_fraction_broadcasted(total, direct)
-        return @. compute_diffuse_fraction(total, direct)
-    end
-
-    frac_diff = TimeVaryingInput(
-        [era5_ncdata_path, era5_ncdata_path],
-        ["msdwswrf", "msdrswrf"],
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        method = time_interpolation_method,
-        compose_function = compute_diffuse_fraction_broadcasted,
-    )
-    LW_d = TimeVaryingInput(
-        era5_ncdata_path,
-        "msdwlwrf",
-        surface_space;
-        start_date,
-        regridder_type,
-        regridder_kwargs = (; interpolation_method),
-        method = time_interpolation_method,
-    )
-    zenith_angle =
-        (t, s) -> default_zenith_angle(
-            t,
-            s;
-            latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
-            longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
-            insol_params = earth_param_set.insol_params,
-        )
-
-    radiation = PrescribedRadiativeFluxes(
-        FT,
-        SW_d,
-        LW_d,
-        start_date;
-        θs = zenith_angle,
-        earth_param_set = earth_param_set,
-        frac_diff = frac_diff,
-    )
-    return (; atmos, radiation)
-end
-
-
 """
      prescribed_lai_era5(era5_lai_ncdata_path,
                          era5_lai_cover_ncdata_path,
diff --git a/src/shared_utilities/perturbed_drivers.jl b/src/shared_utilities/perturbed_drivers.jl
new file mode 100644
index 0000000000..7e4148ecf9
--- /dev/null
+++ b/src/shared_utilities/perturbed_drivers.jl
@@ -0,0 +1,458 @@
+export prescribed_perturbed_temperature_era5, prescribed_perturbed_rh_era5
+
+"""
+    perturbed_temp_specific_humidity_from_dewpoint(dewpoint_temperature, temperature, air_pressure, earth_param_set, ΔT)
+
+Estimates the perturbed specific humidity given the true dewpoint temperature, true temperature of the air
+in Kelvin, and true air pressure in Pa, along with a perturbation to the temperature ΔT, and the ClimaLand
+`earth_param_set`. We first compute relative humidity, and then using RH and the perturbed air temperature,
+equal to  true temperature+ΔT, compute a perturbed `q`.
+
+Relative humidity is computed using the Magnus formula.
+
+For more information on the Magnus formula, see e.g.
+Lawrence, Mark G. (1 February 2005).
+"The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air:
+A Simple Conversion and Applications".
+Bulletin of the American Meteorological Society. 86 (2): 225–234.
+"""
+function perturbed_temp_specific_humidity_from_dewpoint(
+    T_dew_air::data_FT,
+    T_air::data_FT,
+    P_air::data_FT,
+    earth_param_set,
+    ΔT,
+) where {data_FT <: Real}
+    thermo_params = LP.thermodynamic_parameters(earth_param_set)
+    _T_freeze = LP.T_freeze(earth_param_set)
+    sim_FT = typeof(_T_freeze)
+    # Obtain the relative humidity. This function requires temperatures in Celsius
+    rh::sim_FT = rh_from_dewpoint(
+        sim_FT(T_dew_air) - _T_freeze,
+        sim_FT(T_air) - _T_freeze,
+    )
+    q = Thermodynamics.q_vap_from_RH_liquid(
+        thermo_params,
+        sim_FT(P_air),
+        sim_FT(T_air + ΔT),
+        rh,
+    )
+    return q
+end
+
+
+"""
+    perturbed_rh_specific_humidity_from_dewpoint(dewpoint_temperature, temperature, air_pressure, earth_param_set, Δrh)
+
+Estimates the perturbed specific humidity given the true dewpoint temperature, true temperature of the air
+in Kelvin, and true air pressure in Pa, along with a perturbation to the relative humidity Δrh, and the ClimaLand
+`earth_param_set`. We first compute relative humidity, and then perturb it to rh + Δrh, compute a perturbed `q`.
+
+Relative humidity is computed using the Magnus formula.
+
+For more information on the Magnus formula, see e.g.
+Lawrence, Mark G. (1 February 2005).
+"The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air:
+A Simple Conversion and Applications".
+Bulletin of the American Meteorological Society. 86 (2): 225–234.
+"""
+function perturbed_rh_specific_humidity_from_dewpoint(
+    T_dew_air::data_FT,
+    T_air::data_FT,
+    P_air::data_FT,
+    earth_param_set,
+    Δrh,
+) where {data_FT <: Real}
+    thermo_params = LP.thermodynamic_parameters(earth_param_set)
+    _T_freeze = LP.T_freeze(earth_param_set)
+    sim_FT = typeof(_T_freeze)
+    # Obtain the relative humidity. This function requires temperatures in Celsius
+    rh_true::sim_FT = rh_from_dewpoint(
+        sim_FT(T_dew_air) - _T_freeze,
+        sim_FT(T_air) - _T_freeze,
+    )
+    rh = min(max(rh_true + Δrh, sqrt(eps(sim_FT))), sim_FT(1))
+    q = Thermodynamics.q_vap_from_RH_liquid(
+        thermo_params,
+        sim_FT(P_air),
+        sim_FT(T_air),
+        rh,
+    )
+    return q
+end
+
+"""
+     prescribed_perturbed_temperature_era5(era5_ncdata_path,
+                             surface_space,
+                             start_date,
+                             earth_param_set,
+                             ΔT,
+                             FT;
+                             gustiness=1,
+                             max_wind_speed = nothing,
+                             c_co2 = TimeVaryingInput((t) -> 4.2e-4),
+                             time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+                             regridder_type = :InterpolationsRegridder,
+                             interpolation_method = Interpolations.Constant(),)
+
+A helper function which constructs the `PrescribedAtmosphere` and `PrescribedRadiativeFluxes`
+from a file path pointing to the ERA5 data in a netcdf file, the surface_space, the start date,
+and the earth_param_set, applying a change in the instantaneous temperature at each point
+of ΔT, while keeping the relative humidity fixed.
+
+Please see the documentation for `prescribed_forcing_era5` for more information.
+"""
+function prescribed_perturbed_temperature_era5(
+    era5_ncdata_path,
+    surface_space,
+    start_date,
+    earth_param_set,
+    ΔT,
+    FT;
+    gustiness = 1,
+    max_wind_speed = nothing,
+    c_co2 = TimeVaryingInput((t) -> 4.2e-4),
+    time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+    regridder_type = :InterpolationsRegridder,
+    interpolation_method = Interpolations.Constant(),
+)
+    # Pass a list of files in all cases
+    era5_ncdata_path isa String && (era5_ncdata_path = [era5_ncdata_path])
+    _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
+    precip = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["mtpr", "msr"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
+        method = time_interpolation_method,
+        compose_function = (mtpr, msr) -> min.(mtpr .- msr, Float32(0)),
+    )
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
+    snow_precip = TimeVaryingInput(
+        era5_ncdata_path,
+        "msr",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
+        method = time_interpolation_method,
+    )
+    if max_wind_speed isa Nothing
+        wind_function = (u, v) -> sqrt.(u .^ 2 .+ v .^ 2)
+    else
+        wind_function = (u, v) -> min.(sqrt.(u .^ 2 .+ v .^ 2), max_wind_speed)
+    end
+
+    u_atmos = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["u10", "v10"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        compose_function = wind_function,
+        method = time_interpolation_method,
+    )
+    specific_humidity(Td, T, P; params = earth_param_set, ΔT = ΔT) =
+        ClimaLand.perturbed_temp_specific_humidity_from_dewpoint.(
+            Td,
+            T,
+            P,
+            params,
+            ΔT,
+        )
+    q_atmos = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path, era5_ncdata_path],
+        ["d2m", "t2m", "sp"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        compose_function = specific_humidity,
+        method = time_interpolation_method,
+    )
+    P_atmos = TimeVaryingInput(
+        era5_ncdata_path,
+        "sp",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    perturb_temp(T; ΔT = ΔT) = T .+ ΔT
+    T_atmos = TimeVaryingInput(
+        era5_ncdata_path,
+        "t2m",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = perturb_temp),
+        method = time_interpolation_method,
+    )
+    h_atmos = FT(10)
+
+    atmos = PrescribedAtmosphere(
+        precip,
+        snow_precip,
+        T_atmos,
+        u_atmos,
+        q_atmos,
+        P_atmos,
+        start_date,
+        h_atmos,
+        earth_param_set;
+        gustiness = FT(gustiness),
+        c_co2 = c_co2,
+    )
+
+    SW_d = TimeVaryingInput(
+        era5_ncdata_path,
+        "msdwswrf",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    function compute_diffuse_fraction(total, direct)
+        diff = max(total - direct, Float32(0))
+        return min(diff / (total + eps(Float32)), Float32(1))
+    end
+    function compute_diffuse_fraction_broadcasted(total, direct)
+        return @. compute_diffuse_fraction(total, direct)
+    end
+
+    frac_diff = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["msdwswrf", "msdrswrf"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+        compose_function = compute_diffuse_fraction_broadcasted,
+    )
+    LW_d = TimeVaryingInput(
+        era5_ncdata_path,
+        "msdwlwrf",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    zenith_angle =
+        (t, s) -> default_zenith_angle(
+            t,
+            s;
+            latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
+            longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
+            insol_params = earth_param_set.insol_params,
+        )
+
+    radiation = PrescribedRadiativeFluxes(
+        FT,
+        SW_d,
+        LW_d,
+        start_date;
+        θs = zenith_angle,
+        earth_param_set = earth_param_set,
+        frac_diff = frac_diff,
+    )
+    return (; atmos, radiation)
+end
+
+
+"""
+     prescribed_perturbed_rh_era5(era5_ncdata_path,
+                             surface_space,
+                             start_date,
+                             earth_param_set,
+                             Δrh,
+                             FT;
+                             gustiness=1,
+                             max_wind_speed = nothing,
+                             c_co2 = TimeVaryingInput((t) -> 4.2e-4),
+                             time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+                             regridder_type = :InterpolationsRegridder,
+                             interpolation_method = Interpolations.Constant(),)
+
+A helper function which constructs the `PrescribedAtmosphere` and `PrescribedRadiativeFluxes`
+from a file path pointing to the ERA5 data in a netcdf file, the surface_space, the start date,
+and the earth_param_set, applying a change in the instantaneous change to relative
+humidity at each point
+of Δrh. The perturbed rh is clipped to be within the range (0,1].
+
+Please see the documentation for `prescribed_forcing_era5` for more information.
+"""
+function prescribed_perturbed_rh_era5(
+    era5_ncdata_path,
+    surface_space,
+    start_date,
+    earth_param_set,
+    Δrh,
+    FT;
+    gustiness = 1,
+    max_wind_speed = nothing,
+    c_co2 = TimeVaryingInput((t) -> 4.2e-4),
+    time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+    regridder_type = :InterpolationsRegridder,
+    interpolation_method = Interpolations.Constant(),
+)
+    # Pass a list of files in all cases
+    era5_ncdata_path isa String && (era5_ncdata_path = [era5_ncdata_path])
+    _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
+    precip = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["mtpr", "msr"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
+        method = time_interpolation_method,
+        compose_function = (mtpr, msr) -> min.(mtpr .- msr, Float32(0)),
+    )
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
+    snow_precip = TimeVaryingInput(
+        era5_ncdata_path,
+        "msr",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
+        method = time_interpolation_method,
+    )
+    if max_wind_speed isa Nothing
+        wind_function = (u, v) -> sqrt.(u .^ 2 .+ v .^ 2)
+    else
+        wind_function = (u, v) -> min.(sqrt.(u .^ 2 .+ v .^ 2), max_wind_speed)
+    end
+
+    u_atmos = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["u10", "v10"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        compose_function = wind_function,
+        method = time_interpolation_method,
+    )
+    specific_humidity(Td, T, P; params = earth_param_set, Δrh = Δrh) =
+        ClimaLand.perturbed_rh_specific_humidity_from_dewpoint.(
+            Td,
+            T,
+            P,
+            params,
+            Δrh,
+        )
+    q_atmos = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path, era5_ncdata_path],
+        ["d2m", "t2m", "sp"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        compose_function = specific_humidity,
+        method = time_interpolation_method,
+    )
+    P_atmos = TimeVaryingInput(
+        era5_ncdata_path,
+        "sp",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    T_atmos = TimeVaryingInput(
+        era5_ncdata_path,
+        "t2m",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    h_atmos = FT(10)
+
+    atmos = PrescribedAtmosphere(
+        precip,
+        snow_precip,
+        T_atmos,
+        u_atmos,
+        q_atmos,
+        P_atmos,
+        start_date,
+        h_atmos,
+        earth_param_set;
+        gustiness = FT(gustiness),
+        c_co2 = c_co2,
+    )
+
+    SW_d = TimeVaryingInput(
+        era5_ncdata_path,
+        "msdwswrf",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    function compute_diffuse_fraction(total, direct)
+        diff = max(total - direct, Float32(0))
+        return min(diff / (total + eps(Float32)), Float32(1))
+    end
+    function compute_diffuse_fraction_broadcasted(total, direct)
+        return @. compute_diffuse_fraction(total, direct)
+    end
+
+    frac_diff = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["msdwswrf", "msdrswrf"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+        compose_function = compute_diffuse_fraction_broadcasted,
+    )
+    LW_d = TimeVaryingInput(
+        era5_ncdata_path,
+        "msdwlwrf",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    zenith_angle =
+        (t, s) -> default_zenith_angle(
+            t,
+            s;
+            latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
+            longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
+            insol_params = earth_param_set.insol_params,
+        )
+
+    radiation = PrescribedRadiativeFluxes(
+        FT,
+        SW_d,
+        LW_d,
+        start_date;
+        θs = zenith_angle,
+        earth_param_set = earth_param_set,
+        frac_diff = frac_diff,
+    )
+    return (; atmos, radiation)
+end