diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index 3942071af0..60ea0eccb5 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -1647,24 +1647,32 @@ add_diagnostic_variable!(
 function compute_mslp!(out, state, cache, time)
     thermo_params = CAP.thermodynamics_params(cache.params)
     g = TD.Parameters.grav(thermo_params)
-    R_m_surf = Fields.level(
-        lazy.(TD.gas_constant_air.(thermo_params, cache.precomputed.ᶜts)),
-        1,
-    )
+    ts_level = Fields.level(cache.precomputed.ᶜts, 1)
+    R_m_surf = @. lazy(TD.gas_constant_air(thermo_params, ts_level))
 
-    # get pressure, temperature, and height at the lowest atmospheric level
     p_level = Fields.level(cache.precomputed.ᶜp, 1)
-    t_level = Fields.level(
-        lazy.(TD.air_temperature.(thermo_params, cache.precomputed.ᶜts)),
-        1,
-    )
+    t_level = @. lazy(TD.air_temperature(thermo_params, ts_level))
     z_level = Fields.level(Fields.coordinate_field(state.c.ρ).z, 1)
 
-    # compute sea level pressure using the hypsometric equation
+    # Reduce to mean sea level using hypsometric formulation with lapse rate adjustment
+    # Using constant lapse rate Γ = 6.5 K/km, with virtual temperature
+    # represented via R_m_surf. This reduces biases over
+    # very cold or very warm high-topography regions.
+    FT = Spaces.undertype(Fields.axes(state.c.ρ))
+    Γ = FT(6.5e-3) # K m^-1
+
+    #   p_msl = p_z0 * [1 + Γ * z / T_z0]^( g / (R_m Γ))
+    # where:
+    #   - p_z0 pressure at the lowest model level
+    #   - T_z0 air temperature at the lowest model level
+    #   - R_m moist-air gas constant at the surface (R_m_surf), which
+    #     accounts for virtual-temperature effects in the exponent
+    #   - Γ constant lapse rate (6.5 K/km here)
+
     if isnothing(out)
-        return @. p_level * exp(g * z_level / (R_m_surf * t_level))
+        return p_level .* (1 .+ Γ .* z_level ./ t_level) .^ (g / Γ ./ R_m_surf)
     else
-        @. out = p_level * exp(g * z_level / (R_m_surf * t_level))
+        out .= p_level .* (1 .+ Γ .* z_level ./ t_level) .^ (g / Γ ./ R_m_surf)
     end
 end
 
@@ -1673,7 +1681,7 @@ add_diagnostic_variable!(
     long_name = "Mean Sea Level Pressure",
     standard_name = "mean_sea_level_pressure",
     units = "Pa",
-    comments = "Mean sea level pressure computed from the hypsometric equation",
+    comments = "Mean sea level pressure computed using a lapse-rate-dependent hypsometric reduction (ERA-style; Γ=6.5 K/km with virtual temperature via moist gas constant).",
     compute! = compute_mslp!,
 )
 
diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl
index b3708024b5..58be99c97c 100644
--- a/src/initial_conditions/initial_conditions.jl
+++ b/src/initial_conditions/initial_conditions.jl
@@ -386,6 +386,92 @@ function overwrite_initial_conditions!(
     )
 end
 
+"""
+    correct_surface_pressure_for_topography!(
+        p_sfc,
+        file_path,
+        face_space,
+        Y,
+        ᶜT,
+        ᶜq_tot,
+        thermo_params,
+        regridder_kwargs;
+        surface_altitude_var = "z_sfc",
+    )
+
+Adjusts the surface pressure field `p_sfc` to account for mismatches between
+ERA5 (file) surface altitude and the model orography when specifying pressure.
+
+    Δz = z_model_surface - z_sfc
+
+and applies a hydrostatic correction at the surface using the local moist gas
+constant and temperature at the surface:
+
+    p_sfc .= p_sfc .* exp.(-Δz * g ./ (R_m_sfc .* T_sfc))
+
+where:
+- `g` is gravitational acceleration from `thermo_params`
+- `R_m_sfc` is the moist-air gas constant evaluated from `ᶜq_tot` at the surface
+- `T_sfc` is the air temperature from `ᶜT` at the surface
+
+Returns `true` if the correction is applied; returns `false` if the surface
+altitude field cannot be loaded.
+
+Arguments
+- `p_sfc`: face field of surface pressure to be corrected (modified in-place)
+- `file_path`: path to the ERA5-derived initialization NetCDF file
+- `face_space`: face space of the model grid (for reading/regridding)
+- `Y`: prognostic state, used to obtain model surface height
+- `ᶜT`: center field of temperature
+- `ᶜq_tot`: center field of total specific humidity
+- `thermo_params`: thermodynamics parameter set
+- `regridder_kwargs`: keyword arguments forwarded to the regridder
+- `surface_altitude_var`: variable name for surface altitude (default `"z_sfc"`)
+"""
+function correct_surface_pressure_for_topography!(
+    p_sfc,
+    file_path,
+    face_space,
+    Y,
+    ᶜT,
+    ᶜq_tot,
+    thermo_params,
+    regridder_kwargs;
+    surface_altitude_var = "z_sfc",
+)
+    ᶠz_surface = Fields.level(
+        SpaceVaryingInputs.SpaceVaryingInput(
+            file_path,
+            surface_altitude_var,
+            face_space,
+            regridder_kwargs = regridder_kwargs,
+        ),
+        Fields.half,
+    )
+
+    if ᶠz_surface === nothing
+        return false
+    end
+
+    FT = eltype(thermo_params)
+    grav = thermo_params.grav
+
+    ᶠz_model_surface = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+    ᶠΔz = zeros(face_space)
+    @. ᶠΔz = ᶠz_model_surface - ᶠz_surface
+
+    ᶠR_m = ᶠinterp.(TD.gas_constant_air.(thermo_params, TD.PhasePartition.(ᶜq_tot)))
+    ᶠR_m_sfc = Fields.level(ᶠR_m, Fields.half)
+
+    ᶠT = ᶠinterp.(ᶜT)
+    ᶠT_sfc = Fields.level(ᶠT, Fields.half)
+
+    @. p_sfc = p_sfc * exp(FT(-1) * ᶠΔz * grav / (ᶠR_m_sfc * ᶠT_sfc))
+
+    @info "Adjusted surface pressure to account for ERA5/model surface-height differences."
+    return true
+end
+
 # WeatherModel function using the shared implementation
 function overwrite_initial_conditions!(
     initial_condition::WeatherModel,
@@ -436,6 +522,26 @@ function overwrite_initial_conditions!(
         center_space,
         regridder_kwargs = regridder_kwargs,
     )
+    # Apply hydrostatic surface-pressure correction only if surface altitude is available
+    surface_altitude_var = "z_sfc"
+    has_surface_altitude = NC.NCDataset(file_path) do ds
+        haskey(ds, surface_altitude_var)
+    end
+    if has_surface_altitude
+        correct_surface_pressure_for_topography!(
+            p_sfc,
+            file_path,
+            face_space,
+            Y,
+            ᶜT,
+            ᶜq_tot,
+            thermo_params,
+            regridder_kwargs;
+            surface_altitude_var = surface_altitude_var,
+        )
+    else
+        @warn "Skipping topographic correction because variable `$surface_altitude_var` is missing from $(file_path)."
+    end
 
     # With the known temperature (ᶜT) and moisture (ᶜq_tot) profile,
     # recompute the pressure levels assuming hydrostatic balance is maintained.
diff --git a/src/utils/weather_model.jl b/src/utils/weather_model.jl
index b783fa6e91..331a3db0b3 100644
--- a/src/utils/weather_model.jl
+++ b/src/utils/weather_model.jl
@@ -158,13 +158,29 @@ function to_z_levels(source_file, target_file, target_levels, FT)
     end
 
     # Write 2D surface variables - extend to all levels (TODO: accept 2D variables in atmos)
-    # Simply repeat the surface values for all levels
-    surf_map = Dict("skt" => "skt", "sp" => "p")
+    # Duplicate 2D surface field across all target vertical levels
+    surf_map = Dict("skt" => "skt", "sp" => "p", "surface_geopotential" => "z_sfc")
     for (src_name, dst_name) in surf_map
-        var_obj =
-            defVar(ncout, dst_name, FT, ("lon", "lat", "z"), attrib = ncin[src_name].attrib)
+        # Choose attributes; for z_sfc, set clean altitude attributes
+        var_attrib = if dst_name == "z_sfc"
+            Dict(
+                "standard_name" => "surface_altitude",
+                "long_name" => "surface altitude derived from ERA5",
+                "units" => "m",
+                "source_variable" => src_name,
+            )
+        else
+            ncin[src_name].attrib
+        end
+        var_obj = defVar(ncout, dst_name, FT, ("lon", "lat", "z"), attrib = var_attrib)
+        # Read first time slice and coalesce; follow same convention as sp (use [:, :, 1])
+        data2d = FT.(coalesce.(ncin[src_name][:, :, 1], NaN))
+        # Convert geopotential to meters if necessary
+        if dst_name == "z_sfc"
+            data2d .= data2d ./ grav
+        end
         for k in 1:length(target_levels)
-            var_obj[:, :, k] = FT.(ncin[src_name][:, :, 1])
+            var_obj[:, :, k] = data2d
         end
     end