diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 8e6d6bdcfa..c5547111b6 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -549,6 +549,7 @@
 			<var name="soilcl2"/>
 			<var name="soilcl3"/>
 			<var name="soilcl4"/>
+			<var name="soil_color"/>
 #endif
 
                 </stream>
diff --git a/src/core_atmosphere/physics/Registry_noahmp.xml b/src/core_atmosphere/physics/Registry_noahmp.xml
index 89d980f724..b09b57f51c 100644
--- a/src/core_atmosphere/physics/Registry_noahmp.xml
+++ b/src/core_atmosphere/physics/Registry_noahmp.xml
@@ -44,6 +44,10 @@
                 description="soil texture class level 4 needed as input for the NOAH-MP land surface model"
                 packages="sf_noahmp_in"/>
 
+      <var name="soil_color" type="real" dimensions="nCells" units="unitless"
+                description="soil color class"
+                packages="sf_noahmp_in"/>
+
  </var_struct>
 
 
diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F
index 19ad287d00..1f19dedb8e 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F
@@ -235,6 +235,7 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu
 
  integer,pointer:: nsoilcomps
  integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: soil_color
  integer,dimension(:),pointer:: isnowxy
  integer,dimension(:),pointer:: irnumsi,irnummi,irnumfi
  
@@ -323,10 +324,12 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu
  call mpas_pool_get_array(sfc_input,'dzs'   ,dzs   )
  call mpas_pool_get_array(sfc_input,'isltyp',isltyp)
  call mpas_pool_get_array(sfc_input,'ivgtyp',ivgtyp)
+ call mpas_pool_get_array(sfc_input,'soil_color',soil_color)
 
  do i = its, ite
     mpas_noahmp%isltyp(i) = isltyp(i)
     mpas_noahmp%ivgtyp(i) = ivgtyp(i)
+    mpas_noahmp%soil_color(i) = soil_color(i)
  enddo
  do ns = 1, nsoil
     mpas_noahmp%dzs(ns) = dzs(ns,its)
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
index 2de35ed9c2..288d77f098 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
@@ -62,7 +62,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     ! config domain variable
     noahmp%config%domain%SurfaceType                 = 1
     noahmp%config%domain%NumSwRadBand                = 2
-    noahmp%config%domain%SoilColor                   = 4
+    noahmp%config%domain%SoilColor                   = NoahmpIO%soil_color(I)
     noahmp%config%domain%NumCropGrowStage            = 8
     noahmp%config%domain%FlagSoilProcess             = NoahmpIO%calculate_soil
     noahmp%config%domain%NumSoilTimeStep             = NoahmpIO%soil_update_steps
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/EnergyVarOutTransferMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/EnergyVarOutTransferMod.F90
index a3124b3a2e..a150bf9b6e 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/EnergyVarOutTransferMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/EnergyVarOutTransferMod.F90
@@ -136,7 +136,7 @@ subroutine EnergyVarOutTransfer(noahmp, NoahmpIO)
     NoahmpIO%CHB2XY  (I) = noahmp%energy%state%ExchCoeffSh2mBare
     NoahmpIO%Q2MVXY  (I) = noahmp%energy%state%SpecHumidity2mVeg /(1.0-noahmp%energy%state%SpecHumidity2mVeg)  ! spec humidity to mixing ratio
     NoahmpIO%Q2MBXY  (I) = noahmp%energy%state%SpecHumidity2mBare/(1.0-noahmp%energy%state%SpecHumidity2mBare)
-    NoahmpIO%Q2MXY   (I) = noahmp%energy%state%SpecHumidity2m/(1.0-noahmp%energy%state%SpecHumidity2m)
+    NoahmpIO%Q2MXY   (I) = NoahmpIO%Q2MBXY(I) * ( 1 - NoahmpIO%FVEGXY(I) ) + NoahmpIO%Q2MVXY(I) * NoahmpIO%FVEGXY(I)
     NoahmpIO%IRRSPLH (I) = NoahmpIO%IRRSPLH(I) + &
                              (noahmp%energy%flux%HeatLatentIrriEvap * noahmp%config%domain%MainTimeStep)
     NoahmpIO%TSLB    (I,1:NumSoilLayer)       = noahmp%energy%state%TemperatureSoilSnow(1:NumSoilLayer)
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90
index 66458d638e..2782d8bf49 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90
@@ -37,6 +37,7 @@ subroutine NoahmpIOVarFinalizeDefault(NoahmpIO)
     if ( allocated (NoahmpIO%zsoil)     ) deallocate ( NoahmpIO%zsoil               ) ! depth to soil interfaces [m]
     if ( allocated (NoahmpIO%ivgtyp)    ) deallocate ( NoahmpIO%ivgtyp              ) ! vegetation type
     if ( allocated (NoahmpIO%isltyp)    ) deallocate ( NoahmpIO%isltyp              ) ! soil type
+    if ( allocated (NoahmpIO%soil_color)) deallocate ( NoahmpIO%soil_color          ) ! soil color
     if ( allocated (NoahmpIO%vegfra)    ) deallocate ( NoahmpIO%vegfra              ) ! vegetation fraction []
     if ( allocated (NoahmpIO%tmn)       ) deallocate ( NoahmpIO%tmn                 ) ! deep soil temperature [K]
     if ( allocated (NoahmpIO%xland)     ) deallocate ( NoahmpIO%xland               ) ! =2 ocean; =1 land/seaice
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90
index 4f3c3f4f2c..b7fa7077c1 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90
@@ -41,6 +41,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     if ( .not. allocated (NoahmpIO%zsoil)     ) allocate ( NoahmpIO%zsoil      (1:nsoil            )      ) ! depth to soil interfaces [m]
     if ( .not. allocated (NoahmpIO%ivgtyp)    ) allocate ( NoahmpIO%ivgtyp     (its:ite            )      ) ! vegetation type
     if ( .not. allocated (NoahmpIO%isltyp)    ) allocate ( NoahmpIO%isltyp     (its:ite            )      ) ! soil type
+    if ( .not. allocated (NoahmpIO%soil_color)) allocate ( NoahmpIO%soil_color (its:ite            )      ) ! soil color
     if ( .not. allocated (NoahmpIO%vegfra)    ) allocate ( NoahmpIO%vegfra     (its:ite            )      ) ! vegetation fraction []
     if ( .not. allocated (NoahmpIO%tmn)       ) allocate ( NoahmpIO%tmn        (its:ite            )      ) ! deep soil temperature [K]
     if ( .not. allocated (NoahmpIO%xland)     ) allocate ( NoahmpIO%xland      (its:ite            )      ) ! =2 ocean; =1 land/seaice
@@ -469,6 +470,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     NoahmpIO%ice             = undefined_int
     NoahmpIO%ivgtyp          = undefined_int
     NoahmpIO%isltyp          = undefined_int
+    NoahmpIO%soil_color      = undefined_int
     NoahmpIO%isnowxy         = undefined_int
     NoahmpIO%coszen          = undefined_real
     NoahmpIO%xlat            = undefined_real
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90
index 05a29d7035..2b228c571e 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90
@@ -117,6 +117,7 @@ module NoahmpIOVarType
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  soilcl2             ! Soil texture class with depth
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  soilcl3             ! Soil texture class with depth
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  soilcl4             ! Soil texture class with depth
+    real(kind=kind_noahmp), allocatable, dimension(:)      ::  soil_color          ! Soil color class
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  bexp_3D             ! C-H B exponent
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  smcdry_3D           ! Soil Moisture Limit: Dry
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  smcwlt_3D           ! Soil Moisture Limit: Wilt
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index fcbee95775..fda3f4cd79 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -492,6 +492,7 @@
                         <var name="shdmax" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="albedo12m" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="lai12m" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
+                        <var name="soil_color" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soilcomp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soilcl1" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soilcl2" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
@@ -600,6 +601,7 @@
                         <var name="shdmax"/>
                         <var name="albedo12m"/>
                         <var name="lai12m"/>
+                        <var name="soil_color"/>
                         <var name="soilcomp"/>
                         <var name="soilcl1"/>
                         <var name="soilcl2"/>
@@ -931,6 +933,10 @@
                 <var name="iswater_lu" type="integer" dimensions="" units="unitless" default_value="16"
                      description="Index category for water"/>
 
+                <!-- soil color index needed for the Noah-MP land surface scheme -->
+                <var name="soil_color" type="integer" dimensions="nCells" units="unitless"
+                     description="soil color index"/>
+
                 <!-- SOIL COMPOSITION fields needed for the NOAH-MP land surface scheme -->
                 <var name="soilcomp" type="real" dimensions="nSoilComps nCells" units="percent"
                      description="soil composition needed as input in the NOAH-MP land surface model"
diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F
index 2532b4befa..d278050328 100644
--- a/src/core_init_atmosphere/mpas_init_atm_static.F
+++ b/src/core_init_atmosphere/mpas_init_atm_static.F
@@ -1210,6 +1210,16 @@ subroutine init_atm_static(mesh, dims, configs)
 
  if (config_noahmp_static) then
 
+    !
+    ! Interpolate soil_color
+    !
+    geog_sub_path = 'soil_color/'
+
+    call mpas_log_write('--- start interpolate soil_color')
+    call interp_soil_color(mesh, tree, trim(geog_data_path)//trim(geog_sub_path), &
+                          supersample_fac=supersample_fac)
+    call mpas_log_write('--- end interpolate soil_color')
+
     !
     ! Interpolate SOILCOMP
     !
@@ -1777,6 +1787,74 @@ subroutine interp_terrain(mesh, kdtree, geog_data_path, supersample_fac)
 
  end subroutine interp_terrain
 
+ !***********************************************************************
+ !
+ !  routine interp_soil_color
+ !
+ !> \brief   Interpolate soil color
+ !> \author  Michael Barlage
+ !> \date    25 February 2025
+ !> \details
+ !> Interpolate soil color by using the init_atm_map_static_data routine and
+ !> accumulating the pixel values into each cell using categorical_interp_accumulation.
+ !>
+ !> following the method outlined for landuse
+ !
+ !-----------------------------------------------------------------------
+ subroutine interp_soil_color(mesh, kdtree, geog_data_path, supersample_fac)
+
+     implicit none  
+
+     ! Input variables
+     type (mpas_pool_type), intent(inout) :: mesh
+     type (mpas_kd_type), pointer, intent(in) :: kdtree
+     character (len=*), intent(in) :: geog_data_path
+     integer, intent(in), optional :: supersample_fac
+
+     ! Local variables
+     type (mpas_geotile_mgr_type) :: mgr
+     integer, pointer :: nCells
+     integer, dimension(:), pointer :: soil_color
+
+     integer :: iCell
+     integer :: ierr
+
+     ierr = mgr % init(trim(geog_data_path))
+     if (ierr /= 0) then
+         call mpas_log_write("Error occured initalizing interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT)
+         return ! Program execution should not reach this statement since the preceding message is a critical error
+     end if
+
+     call mpas_pool_get_dimension(mesh, 'nCells', nCells)
+     call mpas_pool_get_config(mgr % pool, 'category_min', category_min)
+     call mpas_pool_get_config(mgr % pool, 'category_max', category_max)
+     call mpas_pool_get_array(mesh, 'soil_color', soil_color)
+
+     allocate(ncat(category_min:category_max, nCells))
+     ncat(:,:) = 0
+     
+     call init_atm_map_static_data(mesh, mgr, kdtree, categorical_interp_criteria, categorical_interp_accumulation, &
+                                   supersample_fac=supersample_fac)
+     do iCell = 1, nCells
+         ! Because maxloc returns the location of the maximum value of an array as if the
+         ! starting index of the array is 1, and dataset categories do not necessarily start
+         ! at 1, we need to use category_min to ensure the correct category location is chosen.
+         if (all(ncat(:,iCell) == 0)) then
+           soil_color(iCell) = 11  ! all missing, e.g. water, so choose a neutral value
+         else
+           soil_color(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+         end if
+     end do
+     deallocate(ncat)
+
+     ierr = mgr % finalize()
+     if (ierr /= 0) then
+         call mpas_log_write("Error occured finalizing interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT)
+         return ! Program execution should not reach this statement since the preceding message is a critical error
+     end if
+
+ end subroutine interp_soil_color
+
  !***********************************************************************
  !
  !  routine interp_soilcomp