diff --git a/.github/workflows/frontier/submit.sh b/.github/workflows/frontier/submit.sh
index 61509b0e98..98ce9c4231 100644
--- a/.github/workflows/frontier/submit.sh
+++ b/.github/workflows/frontier/submit.sh
@@ -21,7 +21,7 @@ sbatch <<EOT
 #SBATCH -A CFD154                  # charge account
 #SBATCH -N 1                       # Number of nodes required
 #SBATCH -n 8                       # Number of cores required
-#SBATCH -t 01:00:00                # Duration of the job (Ex: 15 mins)
+#SBATCH -t 01:30:00                # Duration of the job (Ex: 15 mins)
 #SBATCH -o$job_slug.out            # Combined output and error messages file
 #SBATCH -q debug                   # Use debug QOS - only one job per user allowed in queue!
 #SBATCH -W                         # Do not exit until the submitted job terminates.
diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp
index 483de492dd..7a590beea0 100644
--- a/src/common/m_checker_common.fpp
+++ b/src/common/m_checker_common.fpp
@@ -316,8 +316,8 @@ contains
         @:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, &
             "sigma is set but surface_tension is not enabled")
 
-        @:PROHIBIT(surface_tension .and. model_eqns /= 3, &
-            "The surface tension model requires model_eqns=3")
+        @:PROHIBIT(surface_tension .and. (model_eqns /= 3 .and. model_eqns /=2), &
+            "The surface tension model requires model_eqns=3 or model_eqns=2")
 
         @:PROHIBIT(surface_tension .and. num_fluids /= 2, &
             "The surface tension model requires num_fluids=2")
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 329153bfdf..e9d8d0e3e8 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -296,7 +296,7 @@ module m_global_parameters
     !! conditions data to march the solution in the physical computational domain
     !! to the next time-step.
 
-    !$acc declare create(sys_size, buff_size, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, b_size, tensor_size, xi_idx, species_idx, B_idx)
+    !$acc declare create(sys_size, buff_size, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, b_size, tensor_size, xi_idx, species_idx, B_idx, c_idx)
 
     integer :: shear_num !! Number of shear stress components
     integer, dimension(3) :: shear_indices !<
@@ -1212,7 +1212,7 @@ contains
         chemxb = species_idx%beg
         chemxe = species_idx%end
 
-        !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe)
+        !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe, c_idx)
         !$acc update device(b_size, xibeg, xiend, tensor_size)
 
         !$acc update device(species_idx)
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index b6967b913f..0c18d17327 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -135,6 +135,7 @@ contains
 
         real(wp) :: pres_IP, coeff
         real(wp), dimension(3) :: vel_IP, vel_norm_IP
+        real(wp) :: c_IP
         real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP
         real(wp), dimension(nb) :: r_IP, v_IP, pb_IP, mv_IP
         real(wp), dimension(nb*nmom) :: nmom_IP
@@ -170,19 +171,19 @@ contains
             !Interpolate primitive variables at image point associated w/ GP
             if (bubbles_euler .and. .not. qbmm) then
                 call s_interpolate_image_point(q_prim_vf, gp, &
-                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
+                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
                                                r_IP, v_IP, pb_IP, mv_IP)
             else if (qbmm .and. polytropic) then
                 call s_interpolate_image_point(q_prim_vf, gp, &
-                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
+                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
                                                r_IP, v_IP, pb_IP, mv_IP, nmom_IP)
             else if (qbmm .and. .not. polytropic) then
                 call s_interpolate_image_point(q_prim_vf, gp, &
-                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
+                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
                                                r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
             else
                 call s_interpolate_image_point(q_prim_vf, gp, &
-                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP)
+                                               alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
             end if
 
             dyn_pres = 0._wp
@@ -194,6 +195,10 @@ contains
                 q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
             end do
 
+            if (surface_tension) then
+                q_prim_vf(c_idx)%sf(j, k, l) = c_IP
+            end if
+
             if (model_eqns /= 4) then
                 ! If in simulation, use acc mixture subroutines
                 if (elasticity) then
@@ -234,6 +239,11 @@ contains
                 q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
             end do
 
+            ! Set color function
+            if (surface_tension) then
+                q_cons_vf(c_idx)%sf(j, k, l) = c_IP
+            end if
+
             ! Set Energy
             if (bubbles_euler) then
                 q_cons_vf(E_idx)%sf(j, k, l) = (1 - alpha_IP(1))*(gamma*pres_IP + pi_inf + dyn_pres)
@@ -309,6 +319,10 @@ contains
                 q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
             end do
 
+            if (surface_tension) then
+                q_prim_vf(c_idx)%sf(j, k, l) = c_IP
+            end if
+
             call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
                                                             alpha_rho_IP, Re_K, j, k, l)
 
@@ -725,16 +739,18 @@ contains
 
     !> Function that uses the interpolation coefficients and the current state
     !! at the cell centers in order to estimate the state at the image point
-    subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
+    subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
         !$acc routine seq
         type(scalar_field), &
             dimension(sys_size), &
             intent(IN) :: q_prim_vf !< Primitive Variables
+
         real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(INOUT) :: pb, mv
 
         type(ghost_point), intent(IN) :: gp
         real(wp), intent(INOUT) :: pres_IP
         real(wp), dimension(3), intent(INOUT) :: vel_IP
+        real(wp), intent(INOUT) :: c_IP
         real(wp), dimension(num_fluids), intent(INOUT) :: alpha_IP, alpha_rho_IP
         real(wp), optional, dimension(:), intent(INOUT) :: r_IP, v_IP, pb_IP, mv_IP
         real(wp), optional, dimension(:), intent(INOUT) :: nmom_IP
@@ -758,6 +774,8 @@ contains
         pres_IP = 0._wp
         vel_IP = 0._wp
 
+        if (surface_tension) c_IP = 0._wp
+
         if (bubbles_euler) then
             r_IP = 0._wp
             v_IP = 0._wp
@@ -801,6 +819,10 @@ contains
                                       q_prim_vf(advxb + l - 1)%sf(i, j, k)
                     end do
 
+                    if (surface_tension) then
+                        c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)
+                    end if
+
                     if (bubbles_euler .and. .not. qbmm) then
                         !$acc loop seq
                         do l = 1, nb
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index 69a8b4a61b..ab383fb090 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -1650,7 +1650,7 @@ contains
                                     end do
                                 end if
 
-                                ! SURFACE TENSION FLUX. need to check
+                                ! COLOR FUNCTION FLUX
                                 if (surface_tension) then
                                     flux_rs${XYZ}$_vf(j, k, l, c_idx) = &
                                         (xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) + &
@@ -2831,6 +2831,15 @@ contains
                                                 s_P*(xi_R - 1._wp))
                                 end do
 
+                                ! COLOR FUNCTION FLUX
+                                if (surface_tension) then
+                                    flux_rs${XYZ}$_vf(j, k, l, c_idx) = &
+                                        xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) &
+                                        *(vel_L(idx1) + s_M*(xi_L - 1._wp)) &
+                                        + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx) &
+                                        *(vel_R(idx1) + s_P*(xi_R - 1._wp))
+                                end if
+
                                 ! REFERENCE MAP FLUX.
                                 if (hyperelasticity) then
                                     !$acc loop seq