diff --git a/lib/microwaveprop/backtest/features.ex b/lib/microwaveprop/backtest/features.ex
index 2cc0d226..00c26c7a 100644
--- a/lib/microwaveprop/backtest/features.ex
+++ b/lib/microwaveprop/backtest/features.ex
@@ -18,7 +18,11 @@ defmodule Microwaveprop.Backtest.Features do
     a raise on one point kills the whole report.
   """
 
+  import Ecto.Query
+
+  alias Microwaveprop.Repo
   alias Microwaveprop.Weather
+  alias Microwaveprop.Weather.HrrrNativeProfile
 
   @doc """
   Minimum refractivity gradient from the nearest HRRR profile.
@@ -81,4 +85,60 @@ defmodule Microwaveprop.Backtest.Features do
       _ -> nil
     end
   end
+
+  @doc """
+  Surface refractivity from the lowest level of the nearest native
+  hybrid-sigma profile.
+
+  Phase 1 sanity check: this should produce numbers comparable to the
+  `naive_gradient` baseline but computed from native-level data. If
+  the sign/magnitude look wrong, the native ingestion pipeline has a
+  bug.
+
+  Uses the ITU-R P.453-14 formula: N = 77.6*P/T + 3.73e5*e/T²
+  where e (water vapor pressure) is derived from specific humidity.
+  """
+  @spec native_surface_refractivity(float, float, DateTime.t()) :: float | nil
+  def native_surface_refractivity(lat, lon, valid_time) do
+    with %HrrrNativeProfile{} = profile <- find_nearest_native(lat, lon, valid_time),
+         t when is_float(t) <- profile.surface_temp_k,
+         p when is_float(p) <- profile.surface_pressure_pa,
+         q when is_float(q) <- profile.surface_spfh do
+      # Water vapor pressure from specific humidity: e = q*P / (0.622 + 0.378*q)
+      e = q * p / (0.622 + 0.378 * q)
+      # N-units
+      77.6 * p / (t * 100) + 3.73e5 * e / (t * t * 100)
+    else
+      _ -> nil
+    end
+  end
+
+  defp find_nearest_native(lat, lon, valid_time) do
+    dlat = 0.07
+    dlon = 0.07
+    time_start = DateTime.add(valid_time, -3600, :second)
+    time_end = DateTime.add(valid_time, 3600, :second)
+
+    HrrrNativeProfile
+    |> where(
+      [p],
+      p.lat >= ^(lat - dlat) and p.lat <= ^(lat + dlat) and
+        p.lon >= ^(lon - dlon) and p.lon <= ^(lon + dlon) and
+        p.valid_time >= ^time_start and p.valid_time <= ^time_end
+    )
+    |> order_by([p],
+      asc:
+        fragment(
+          "ABS(? - ?) + ABS(? - ?) + ABS(EXTRACT(EPOCH FROM ? - ?))",
+          p.lat,
+          ^lat,
+          p.lon,
+          ^lon,
+          p.valid_time,
+          ^valid_time
+        )
+    )
+    |> limit(1)
+    |> Repo.one()
+  end
 end
diff --git a/lib/microwaveprop/weather/grib2/wgrib2.ex b/lib/microwaveprop/weather/grib2/wgrib2.ex
index efa2aadd..cc56af07 100644
--- a/lib/microwaveprop/weather/grib2/wgrib2.ex
+++ b/lib/microwaveprop/weather/grib2/wgrib2.ex
@@ -251,15 +251,23 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do
     points_per_message = nx * ny
     bytes_per_message = points_per_message * 4
 
+    # wgrib2 -lola ... bin writes Fortran unformatted binary: each
+    # message is preceded by a 4-byte little-endian record length
+    # and followed by a duplicate 4-byte record length. So the
+    # stride per message is 4 + data + 4 = data + 8.
+    record_overhead = 8
+    stride = bytes_per_message + record_overhead
+
     result =
       messages
       |> Enum.with_index()
       |> Enum.reduce(%{}, fn {msg, msg_idx}, acc ->
-        offset = msg_idx * bytes_per_message
+        # Skip the 4-byte header to reach the data
+        data_offset = msg_idx * stride + 4
         key = "#{msg.var}:#{msg.level}"
 
-        if offset + bytes_per_message <= byte_size(bin_data) do
-          chunk = binary_part(bin_data, offset, bytes_per_message)
+        if data_offset + bytes_per_message <= byte_size(bin_data) do
+          chunk = binary_part(bin_data, data_offset, bytes_per_message)
           merge_message_values(acc, key, chunk, nx, ny, lon_start, lon_step, lat_start, lat_step)
         else
           acc
diff --git a/lib/microwaveprop/weather/hrrr_native_client.ex b/lib/microwaveprop/weather/hrrr_native_client.ex
index b5fc31a4..faacb0c0 100644
--- a/lib/microwaveprop/weather/hrrr_native_client.ex
+++ b/lib/microwaveprop/weather/hrrr_native_client.ex
@@ -139,7 +139,9 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do
   def extract_native_profiles(grib_binary, points) when is_list(points) do
     alias Microwaveprop.Weather.Grib2.Wgrib2
 
-    match_pattern = ":(#{Enum.join(@native_variables, "|")}):"
+    # Match only hybrid-level messages — the simple var-name pattern
+    # also hits surface/2m/10m messages which misalign the binary output.
+    match_pattern = ":(#{Enum.join(@native_variables, "|")}):.*hybrid level:"
 
     if Wgrib2.available?() do
       grid_spec = bounding_grid(points)
diff --git a/test/microwaveprop/backtest/features_test.exs b/test/microwaveprop/backtest/features_test.exs
index a67d750f..b7430273 100644
--- a/test/microwaveprop/backtest/features_test.exs
+++ b/test/microwaveprop/backtest/features_test.exs
@@ -3,6 +3,7 @@ defmodule Microwaveprop.Backtest.FeaturesTest do
 
   alias Microwaveprop.Backtest.Features
   alias Microwaveprop.Weather
+  alias Microwaveprop.Weather.HrrrNativeProfile
 
   @point_lat 32.9
   @point_lon -97.0
@@ -59,4 +60,37 @@ defmodule Microwaveprop.Backtest.FeaturesTest do
       assert Features.pressure(@point_lat, @point_lon, @valid_time) == nil
     end
   end
+
+  describe "native_surface_refractivity/3" do
+    test "computes N-units from native profile surface scalars" do
+      {:ok, _} =
+        %HrrrNativeProfile{}
+        |> HrrrNativeProfile.changeset(%{
+          valid_time: @valid_time,
+          lat: @point_lat,
+          lon: @point_lon,
+          level_count: 1,
+          heights_m: [10.0],
+          temp_k: [295.0],
+          spfh: [0.010],
+          pressure_pa: [101_325.0],
+          u_wind_ms: [2.0],
+          v_wind_ms: [1.0],
+          tke_m2s2: [0.5],
+          surface_temp_k: 295.0,
+          surface_spfh: 0.010,
+          surface_pressure_pa: 101_325.0
+        })
+        |> Repo.insert()
+
+      result = Features.native_surface_refractivity(@point_lat, @point_lon, @valid_time)
+      assert is_float(result)
+      # Typical surface N is 250-400 in temperate latitudes
+      assert result > 200.0 and result < 500.0
+    end
+
+    test "returns nil when no native profile exists" do
+      assert Features.native_surface_refractivity(@point_lat, @point_lon, @valid_time) == nil
+    end
+  end
 end