diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
index 460b9ecda428..debc65395c74 100755
--- a/components/mpas-ocean/bld/build-namelist
+++ b/components/mpas-ocean/bld/build-namelist
@@ -783,6 +783,8 @@ add_default($nl, 'config_frazil_in_open_ocean');
 add_default($nl, 'config_frazil_under_land_ice');
 add_default($nl, 'config_frazil_heat_of_fusion');
 add_default($nl, 'config_frazil_ice_density');
+add_default($nl, 'config_use_frazil_ice_porosity');
+add_default($nl, 'config_frazil_ice_porosity');
 add_default($nl, 'config_frazil_fractional_thickness_limit');
 add_default($nl, 'config_specific_heat_sea_water');
 add_default($nl, 'config_frazil_maximum_depth');
diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
index c5dad5d935a9..1d09fdfbdb45 100644
--- a/components/mpas-ocean/bld/build-namelist-section
+++ b/components/mpas-ocean/bld/build-namelist-section
@@ -291,6 +291,8 @@ add_default($nl, 'config_frazil_in_open_ocean');
 add_default($nl, 'config_frazil_under_land_ice');
 add_default($nl, 'config_frazil_heat_of_fusion');
 add_default($nl, 'config_frazil_ice_density');
+add_default($nl, 'config_use_frazil_ice_porosity');
+add_default($nl, 'config_frazil_ice_porosity');
 add_default($nl, 'config_frazil_fractional_thickness_limit');
 add_default($nl, 'config_specific_heat_sea_water');
 add_default($nl, 'config_frazil_maximum_depth');
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
index 4f9b5d8fe4d3..559bc1ed84e3 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -394,6 +394,8 @@
 <config_frazil_under_land_ice>.true.</config_frazil_under_land_ice>
 <config_frazil_heat_of_fusion>3.337e5</config_frazil_heat_of_fusion>
 <config_frazil_ice_density>1000.0</config_frazil_ice_density>
+<config_use_frazil_ice_porosity>.false.</config_use_frazil_ice_porosity>
+<config_frazil_ice_porosity>1.0</config_frazil_ice_porosity>
 <config_frazil_fractional_thickness_limit>0.1</config_frazil_fractional_thickness_limit>
 <config_specific_heat_sea_water>3.996e3</config_specific_heat_sea_water>
 <config_frazil_maximum_depth>100.0</config_frazil_maximum_depth>
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
index bea1e98d9de8..09ef82120433 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -1539,6 +1539,22 @@ Valid values: Any positive real number.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_use_frazil_ice_porosity" type="logical"
+	category="frazil_ice" group="frazil_ice">
+If true, use an internal porosity to set internal frazil salinity and that of specified passive tracers.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_frazil_ice_porosity" type="real"
+	category="frazil_ice" group="frazil_ice">
+If frazil porosity is used, sets the porosity of internal ocean frazil ice.
+
+Valid values: Positive real number from 0 to 1.
+Default: Defined in namelist_defaults.xml
+</entry>
+
 <entry id="config_frazil_fractional_thickness_limit" type="real"
 	category="frazil_ice" group="frazil_ice">
 maximum fraction of layer thickness than can be used or created at an instant by frazil.
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 2a21e22f40d2..e4276a705e5d 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -985,10 +985,18 @@
 					description="Energy per kilogram released when sea water freezes. NOTE: test and make consistent with E3SM."
 					possible_values="Any positive real number."
 		/>
+		<nml_option name="config_use_frazil_ice_porosity" type="logical" default_value=".false."
+					description="If true, use an internal porosity to set internal frazil salinity and that of specified passive tracers."
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_frazil_ice_density" type="real" default_value="1000.0" units="kg m^-3"
 					description="Assumed density of frazil. NOTE: test and make consistent with E3SM."
 					possible_values="Any positive real number."
 		/>
+		<nml_option name="config_frazil_ice_porosity" type="real" default_value="1.0" units="1"
+					description="If frazil porosity is used, sets the porosity of internal ocean frazil ice."
+					possible_values="Positive real number from 0 to 1"
+		/>
 		<nml_option name="config_frazil_fractional_thickness_limit" type="real" default_value="0.1" units="non-dimensional"
 					description="maximum fraction of layer thickness than can be used or created at an instant by frazil."
 					possible_values="Any positive real number between 0 and 1."
@@ -4076,6 +4084,10 @@
 			 description="surface pressure forcing due to weight of frazil ice"
 			 packages="frazilIce"
 		/>
+		<var name="lowSalinityFrazilIce" type="real" dimensions="nCells Time" units="kg m^-2"
+			 description="min(0,AccumulatedFrazilIceMass * reference salinity - accumulatedFrazilIceSalinity), non-zero indicates possible non-conservation"
+			 packages="frazilIce"
+		/>
 		<var name="frazilIceFreshwaterFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
 			description="Flux of frazil ice mass through the ocean surface. Negative is defined as mass leaving the ocean."
 			packages="frazilIce"
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
index 5ca55901e3eb..c5d8b7185fe7 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
@@ -360,6 +360,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       real (kind=RKIND), dimension(:,:), pointer :: frazilTemperatureTendency
       real (kind=RKIND), dimension(:,:), pointer :: frazilSalinityTendency
       real (kind=RKIND), dimension(:), pointer :: frazilSurfacePressure
+      real (kind=RKIND), dimension(:), pointer :: lowSalinityFrazilIce
       integer, dimension(:), pointer :: landIceMask, landIceFloatingMask
 
       integer :: iCell, k, nCells
@@ -390,7 +391,9 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       integer, pointer :: indexSalinityPtr !< index in tracers array for salinity
       integer :: indexTemperature !< index in tracers array for temperature
       integer :: indexSalinity !< index in tracers array for salinity
-      integer :: kBottomFrazil  ! k index where testing for frazil begins
+      integer :: kBottomFrazil ! k index where testing for frazil begins
+      integer :: usePorosity   ! 1 when computing internal frazil salinity using porosity
+                               ! 0 when using config_frazil_sea_ice_reference_salinity
 
       logical :: underFloatingLandIce, underLandIce ! indicates if we are under land ice
 
@@ -435,6 +438,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       call mpas_pool_get_array(forcingPool, 'frazilSalinityTendency', frazilSalinityTendency)
       call mpas_pool_get_array(forcingPool, 'frazilLayerThicknessTendency', frazilLayerThicknessTendency)
       call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure)
+      call mpas_pool_get_array(forcingPool, 'lowSalinityFrazilIce', lowSalinityFrazilIce)
       call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask)
       call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask)
       call mpas_pool_get_array(forcingPool, 'frazilIceFreshwaterFlux', frazilIceFreshwaterFlux)
@@ -448,6 +452,13 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       frazilSalinityTendency = 0.0_RKIND
       frazilLayerThicknessTendency = 0.0_RKIND
 
+      ! initialize frazil salinity conservation field
+      lowSalinityFrazilIce = 0.0_RKIND
+
+      ! Determine if using porosity to compute internal frazilSalinity
+      usePorosity = 0 ! use config_frazil_sea_ice_reference_salinity
+      if (config_use_frazil_ice_porosity) usePorosity = 1
+
       ! frazil fields are needed only over 0 and 1 halos
       nCells = nCellsArray( 2 )
 
@@ -489,7 +500,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
         ! find deepest level where frazil can be created
         kBottomFrazil=maxLevelCell(iCell)
         do k=maxLevelCell(iCell), minLevelCell(iCell), -1
-          ! add the ssh so frazil can form below land ice, where the ssh is depressed
+         ! add the ssh so frazil can form below land ice, where the ssh is depressed
           if (-zMid(k,iCell) < -ssh(iCell) + config_frazil_maximum_depth) then
             kBottomFrazil=k
             exit
@@ -529,9 +540,10 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
 
             ! get frazil salinity
             if (underLandIce) then
-              frazilSalinity = 0.0_RKIND
+              frazilSalinity = usePorosity * config_frazil_ice_porosity * activeTracers(indexSalinity, k, iCell)
             else
-              frazilSalinity = config_frazil_sea_ice_reference_salinity
+              frazilSalinity = usePorosity * config_frazil_ice_porosity * activeTracers(indexSalinity, k, iCell) &
+                 + config_frazil_sea_ice_reference_salinity * (1 - usePorosity)
             end if
             frazilSalinity = min( frazilSalinity, activeTracers(indexSalinity, k, iCell) )
 
@@ -541,7 +553,6 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
             ! limit the frazil formed appropriately
             newFrazilIceThickness = min(newFrazilIceThickness, layerThickness(k,iCell) * config_frazil_fractional_thickness_limit)
 
-
             ! Determine new salt in frazil
             newThicknessWeightedSaltContent = newFrazilIceThickness * frazilSalinity
 
@@ -564,7 +575,6 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
             sumNewThicknessWeightedSaltContent = sumNewThicknessWeightedSaltContent + newThicknessWeightedSaltContent
 
           else
-
             ! ocean water is warm enough to melt frazil
 
             ! test to see if there is frazil to be melted
@@ -600,8 +610,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
                                                  / dt
 
               ! keep track of new frazil ice
-              sumNewThicknessWeightedSaltContent = sumNewThicknessWeightedSaltContent - meltedThicknessWeightedSaltContent
-              sumNewFrazilIceThickness = sumNewFrazilIceThickness - meltedFrazilIceThickness
+              sumNewThicknessWeightedSaltContent = max(0.0_RKIND,sumNewThicknessWeightedSaltContent - meltedThicknessWeightedSaltContent)
+              sumNewFrazilIceThickness = max(0.0_RKIND,sumNewFrazilIceThickness - meltedFrazilIceThickness)
 
             endif  ! if (sumNewFrazilIceThickness > 0.0_RKIND)
 
@@ -609,22 +619,42 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
 
         enddo   ! do k=kBottom,minLevelCell,-1
 
-        ! frazilIceFreshwaterFlux forms part of the landIceFreshwaterFluxTotal computed in surface land ice fluxes
-        frazilIceFreshwaterFlux(iCell) = - (sumNewFrazilIceThickness * config_frazil_ice_density) / dt
         ! accumulate frazil mass to column total
         ! note: the accumulatedFrazilIceMass (at both time levels) is reset to zero after being sent to the coupler
         accumulatedFrazilIceMassNew(iCell) = accumulatedFrazilIceMassOld(iCell) + sumNewFrazilIceThickness &
-                                           * config_frazil_ice_density
-        accumulatedFrazilIceSalinityNew(iCell) = accumulatedFrazilIceSalinityOld(iCell) + sumNewThicknessWeightedSaltContent
+           * config_frazil_ice_density
+
+        accumulatedFrazilIceSalinityNew(iCell) = accumulatedFrazilIceSalinityOld(iCell) + &
+           (1-usePorosity) * sumNewThicknessWeightedSaltContent
 
+        !NJ: Assume frazil salinity is the coupling salinity
         if ( underLandIce ) then
+           ! frazilIceFreshwaterFlux forms part of the landIceFreshwaterFluxTotal computed in surface land ice fluxes
+           frazilIceFreshwaterFlux(iCell) = - (sumNewFrazilIceThickness * config_frazil_ice_density) / dt
            ! accumulate frazil formed under land ice in case we're not coupling and we need to keep track of it
            ! for freshwater budgets
            accumulatedLandIceFrazilMassNew(iCell) = accumulatedLandIceFrazilMassOld(iCell) &
                                                   + sumNewFrazilIceThickness * config_frazil_ice_density
            ! There is no accumulatedLandIceFrazilSalinity because it is always 0
+           frazilSalinity = 0.0_RKIND
+        else
+           frazilSalinity = config_frazil_sea_ice_reference_salinity
         end if
 
+        !Determine salt expelled in the upper layer to arrive at the assumed frazil salt for coupling
+        if (sumNewFrazilIceThickness > 0.0_RKIND .and. usePorosity .ge. 1) then
+           newThicknessWeightedSaltContent = sumNewFrazilIceThickness * frazilSalinity
+           if (newThicknessWeightedSaltContent > sumNewThicknessWeightedSaltContent) then
+            !Warning: salt conservation error. Track with lowSalinityFrazilIce
+            lowSalinityFrazilIce(iCell) = newThicknessWeightedSaltContent - sumNewThicknessWeightedSaltContent
+            accumulatedFrazilIceSalinityNew(iCell) = accumulatedFrazilIceSalinityOld(iCell) + max(0.0_RKIND, sumNewThicknessWeightedSaltContent)
+           else
+            !salt expelled at the surface
+            frazilSalinityTendency(minLevelCell(iCell),iCell) = frazilSalinityTendency(minLevelCell(iCell),iCell) + &
+                max(0.0_RKIND,(sumNewThicknessWeightedSaltContent -  newThicknessWeightedSaltContent) ) / dt
+            accumulatedFrazilIceSalinityNew(iCell) = accumulatedFrazilIceSalinityOld(iCell) + newThicknessWeightedSaltContent
+           endif
+         endif ! if frazil ice is formed
       enddo   ! do iCell = 1, nCells
       !$omp end do
       !$omp end parallel