diff --git a/large_media b/large_media
index 8cf736f69045..fb009ce0b606 160000
--- a/large_media
+++ b/large_media
@@ -1 +1 @@
-Subproject commit 8cf736f690459f9829bbda58da5f4a76baa548c5
+Subproject commit fb009ce0b606195537fb3bdb03e338c5c49f39b0
diff --git a/modules/combined/doc/content/modules/combined/examples/stm_laserwelding_dimred.md b/modules/combined/doc/content/modules/combined/examples/stm_laserwelding_dimred.md
new file mode 100644
index 000000000000..159a96c1879e
--- /dev/null
+++ b/modules/combined/doc/content/modules/combined/examples/stm_laserwelding_dimred.md
@@ -0,0 +1,172 @@
+# Using Stochastic Tools for the full-field reconstruction of multiphysics problems
+
+The purpose of this example is to present a multiphysics model order reduction case using
+the [Stochastic Tools Module](modules/stochastic_tools/index.md). The problem of interest is a
+melt pool model using the [Navier-Stokes](modules/navier_stokes/index.md) module.
+The problem involves multiple uncertain process parameters and are interested in predicting
+full field quantities for unseen process parameter combinations using non-intrusive surrogate
+models in combination with Proper Orthogonal Decomposition (POD).
+
+## Problem Description
+
+The problem of interest involves a transient 2D melt pool simulation for a 316L stainless steel
+with a moving Gaussian laser beam. The laser source melts the top layer of the metal and
+based on the equlibrium of the forces (surface tension, vapor recoil pressure) the surface of the
+melt pool can deform considerably. The changing geometry is taken into account using an
+Arbitrary Legrangian Eulerian (ALE) approach. The setup of the problem is depicted in [fig:geom].
+
+!media combined/laser_weld/geometry.png caption=Problem setup for the laser melt pool simulation id=fig:geom style=width:100%
+
+The temperature-dependent thermal properties of the problem have been taken from the
+references listed in [tab:mat_prop].
+
+!table caption=Material properties for 316L Stailess steel id=tab:mat_prop
+| Property | Source|
+| :- | :- |
+| Dynamic viscosity | [!cite](kim1975thermophysical)$^\dagger$ |
+| Density | [!cite](kim1975thermophysical) |
+| Specific heat | [!cite](kim1975thermophysical) |
+| Thermal conductivity | [!cite](pichler2022surface) |
+| Surface tension | [!cite](pichler2022surface) |
+| Vapor recoil pressure | [!cite](chen2021numerical) |
+
+$\dagger$ In the treatment of dynamic viscosity, following the recommendations in [!cite](noble2007use),
+we artificially cap the minimum value to stabilize the fluid flow solver.
+
+The input file used for the simulations is presented below:
+
+!listing examples/stochastic/laser_welding_dimred/2d.i caption=Laser meltpool model input file id=list:meltpool-ref
+
+We see that it utilizes the `!include` syntax to read the process parameters togheter with the mesh and
+kernel and boundary condition objects from other input files. This input file design reduces
+duplication when it comes to the testing of the reduced-order models in subsequent sections.
+
+!listing examples/stochastic/laser_welding_dimred/parameters.i caption=Input file for specifying process parameters id=list:parameters-ref
+
+!listing examples/stochastic/laser_welding_dimred/mesh.i caption=Input file for specifying the mesh id=list:mesh-ref
+
+!listing examples/stochastic/laser_welding_dimred/physics_objects.i caption=Input file adding the objects that define the physics id=list:physics-ref
+
+### Process parameters
+
+A total of two process parameters have been identified to be changable: the effective laser power and
+the effective laser radius. They can change in the intervals shown in [tab:ppars] assuming a uniform
+distribution.
+
+!table caption=The distributions used for the process parameters id=tab:ppars
+| Property | Min. | Max. |
+| :- | - | - |
+| Effective laser power  | 60 W | 74 W |
+| Effective laser radius | 125 $\mu m$ | 155 $\mu m$ |
+
+### Quantities of Interest
+
+The quantity of interest in our case is the whole temperature field at the end of the simulation (0.4~ms).
+To be able to reconstruct the whole field we utilize Proper Orthogonal Decomposition (POD) in conjunction
+ with a Multi-output Gaussian Process
+(MOGP).
+
+## Results
+
+## Reference results
+
+Two reference results are presented here:
+1. Results for a simulation with low effective laser power ($60~W$) and wide effective laser radius ($155~\mu m$)
+2. Results for a simulation with high effective laser power ($74~W$) and narrow effective laser radius ($125~\mu m$)
+
+The reference results are filtered to show the melt pool shape. We see that for the first case, the
+steel melts, but no displacement is visible due to little to no evaporation. On the other hand, the
+increase in laser power results in an increased evaporation resulting in significant surface deformation.
+
+!media combined/laser_weld/60.png caption=Simulation results at the final time step with laser power ($60~W$) and wide effective laser radius ($155~\mu m$) id=fig:result-60 style=width:100%
+
+!media combined/laser_weld/74.png caption=Simulation results at the final time step with laser power ($74~W$) and narrow effective laser radius ($125~\mu m$) id=fig:result-74 style=width:100%
+
+## Training the reduced-order model
+
+The reduced-order model was trained by first running the reference model with 45 different
+process parameter combinations in a [Monte Carlo sampler](MonteCarloSampler.md). The input file used for training
+is presented below.
+
+!listing examples/stochastic/laser_welding_dimred/train.i id=list:training
+
+The parts that require special attention are the snapshot gathering
+parts in the transfer and reporter blocks.
+
+!listing examples/stochastic/laser_welding_dimred/train.i caption=Transfers and Reporters needed for snapshot collection block=Transfers Reporters VariableMappings remove=Reporters/matrix Reporters/svd id=list:training-trans-rep
+
+With the snapshots all available in one place, we can extract the POD modes
+necessary for the dimensionality reduction. This is done using the
+[PODMapping.md] reporter which not only creates the decomposition but also maps the
+created snapshots onto a low-dimensional latent space.
+
+!plot scatter
+  id=results caption=Scree plot of the first 8 singular values extracted from the training data.
+  filename=examples/stochastic/laser_welding_dimred/gold/sv_to_print.csv
+  data=[{'x':'i', 'y':'sv', 'name':'Singular values'}]
+  layout={'xaxis':{'type':'linear', 'title':'Mode #'},
+          'yaxis':{'type':'log','title':'Singular Value (-)'}}
+
+ In this specific example, the
+dimensionality of the latent space was selected to be 8.
+Once the coordinates of the snapshots are available in the latent space, we can fit a
+[MOGP](GaussianProcessTrainer.md) on them. This is done in the `Trainers` block with the corresponding covariance functions defined in the `Covariance` block.
+
+!listing examples/stochastic/laser_welding_dimred/train.i block=Trainers Covariance id=list:training-gp caption=Definition of the Gaussian Process surrogate used for fitting the coordinates in the latent space.
+
+One can then save both the POD mapping and MOGP into restratable structures for
+later usage.
+
+!listing examples/stochastic/laser_welding_dimred/train.i block=Outputs
+
+## Evaluating the error of the ROM
+
+We evaluate the ROM by comparing reconstructed temperature fields at unseen process
+parameter combinations to full-order model results at the same samples.
+The metric for comparison is the relative $L^2$ error defined as:
+
+\begin{equation}
+e = \frac{||T_{ROM}-T_{FOM}||}{||T_{FOM}||}
+\end{equation}
+
+This error is generated in a new input file, which computes the reconstructed field and
+solves the full-order model at the same time.
+
+!listing examples/stochastic/laser_welding_dimred/2d-reconst.i caption=Test input file for comparing the reconstructed fields to the reference solutions id=list:2d-reconst
+
+We see that duplication is minimized by reusing the same building blocks as the original input file.
+In this case, we load both the MOGP and the POD mapping from a file:
+
+!listing examples/stochastic/laser_welding_dimred/2d-reconst.i block=VariableMappings Surrogates
+
+And use them in an object that reconstructs the field variable and inserts it into a `T_reconst`
+auxiliary variable.
+
+!listing examples/stochastic/laser_welding_dimred/2d-reconst.i block=UserObjects
+
+Then, the integrated L2 error is computed by the appropriate postprocessor.
+
+!listing examples/stochastic/laser_welding_dimred/2d-reconst.i block=Postprocessors
+
+This input file is run for 90 samples of new process parameter combination selected by
+a Monte Carlo sampler. The input file for automating the process of error analysis is presented below.
+
+!listing examples/stochastic/laser_welding_dimred/test.i id=list:2d-reconst-pp caption=Input file for the automation of the testing phase.
+
+This simply gathers the relative $L^2$ errors from the subapplications and
+presents the results in a `json` file.
+The error histogram is presented below. We see that the maximum integrated $L^2$
+error is below 1.6% with the majority of the samples having error
+in the 0-0.2% interval.
+
+
+!plot histogram
+  id=test-results caption=Histogram of the L2 errors over the 90 test samples.
+  filename=examples/stochastic/laser_welding_dimred/gold/error_to_plot.csv
+  data=[{'x':'l2error', 'name':'L2 Error'}]
+  layout={'xaxis':{'type':'linear', 'title':'Relative L2 Error'},
+          'yaxis':{'type':'linear','title':'Frequency'}}
+
+
+
+
diff --git a/modules/combined/doc/content/modules/combined/examples/stm_thermomechanics.md b/modules/combined/doc/content/modules/combined/examples/stm_thermomechanics.md
index cf0ddbc8250e..c34deafc5310 100644
--- a/modules/combined/doc/content/modules/combined/examples/stm_thermomechanics.md
+++ b/modules/combined/doc/content/modules/combined/examples/stm_thermomechanics.md
@@ -44,7 +44,7 @@ The problem of interest involves a steady-state thermomechanics model. The geome
 
 !row-end!
 
-!listing examples/stochastic/graphite_ring_thermomechanics.i caption=Thermomechanics model input file id=list:thermo
+!listing examples/stochastic/thermomech/graphite_ring_thermomechanics.i caption=Thermomechanics model input file id=list:thermo
 
 ### Uncertain Parameters
 
@@ -96,11 +96,11 @@ Using [latin hypercube sampling](LatinHypercubeSampler.md), the thermomechanics
 | Polynomial Chaos --- Training | 7,344 | 13.7 hr  |
 | Polynomial Chaos --- Evaluation | 100,000 | 6.8 s  |
 
-!listing combined/examples/stochastic/lhs_uniform.i id=list:lhs caption=Latin hypercube sampling and statistics input file
+!listing combined/examples/stochastic/thermomech/lhs_uniform.i id=list:lhs caption=Latin hypercube sampling and statistics input file
 
-!listing combined/examples/stochastic/poly_chaos_train_uniform.i id=list:train caption=Polynomial chaos training input file
+!listing combined/examples/stochastic/thermomech/poly_chaos_train_uniform.i id=list:train caption=Polynomial chaos training input file
 
-!listing combined/examples/stochastic/poly_chaos_uniform.i id=list:eval caption=Polynomial chaos evaluation input file
+!listing combined/examples/stochastic/thermomech/poly_chaos_uniform.i id=list:eval caption=Polynomial chaos evaluation input file
 
 ### Statistics
 
diff --git a/modules/combined/doc/content/modules/combined/tutorials/index.md b/modules/combined/doc/content/modules/combined/tutorials/index.md
index c970c6f033e2..5075ec991257 100644
--- a/modules/combined/doc/content/modules/combined/tutorials/index.md
+++ b/modules/combined/doc/content/modules/combined/tutorials/index.md
@@ -5,6 +5,7 @@ physics modules can be used together:
 
 - [Introduction to Coupled Thermo-Mechanical Modeling](combined/tutorials/introduction/index.md)
 - [Using Stochastic Tools with Multiphysics Models](combined/examples/stm_thermomechanics.md)
+- [Using dimensionality reduction and full-field reconstruction with multiphysics models](combined/examples/stm_laserwelding_dimred.md)
 
 # Solid Isotropic Material Penalization (SIMP) Topology Optimization
 
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/2d-reconst.i b/modules/combined/examples/stochastic/laser_welding_dimred/2d-reconst.i
new file mode 100644
index 000000000000..7f55905ffac4
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/2d-reconst.i
@@ -0,0 +1,87 @@
+!include parameters.i
+
+!include mesh.i
+
+[Variables]
+  [vel]
+    family = LAGRANGE_VEC
+  []
+  [T]
+  []
+  [p]
+  []
+  [disp_x]
+  []
+  [disp_y]
+  []
+[]
+
+[AuxVariables]
+  [T_reconst]
+  []
+[]
+
+!include physics_objects.i
+
+[UserObjects]
+  [inverse]
+    type = InverseMapping
+    mapping = pod
+    variable_to_fill = "T_reconst"
+    variable_to_reconstruct = "T"
+    surrogate = mogp
+    parameters = '${R} ${power}'
+    execute_on = TIMESTEP_END
+  []
+[]
+
+[Surrogates]
+  [mogp]
+    type = GaussianProcessSurrogate
+    filename = "train_mogp_out_mogp.rd"
+  []
+[]
+
+[VariableMappings]
+  [pod]
+    type = PODMapping
+    filename = "train_pod_out_pod.rd"
+  []
+[]
+
+[Executioner]
+  type = Transient
+  end_time = ${endtime}
+  dtmin = 1e-10
+  dtmax = 1e-5
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu       NONZERO'
+  petsc_options = '-snes_converged_reason -ksp_converged_reason -options_left'
+  solve_type = 'NEWTON'
+  line_search = 'none'
+  nl_max_its = 16
+  l_max_its = 100
+  [TimeStepper]
+    type = IterationAdaptiveDT
+    optimal_iterations = 10
+    iteration_window = 4
+    dt = ${timestep}
+    linear_iteration_ratio = 1e6
+    growth_factor = 1.25
+  []
+[]
+
+[Postprocessors]
+  [l2error]
+    type = ElementL2Difference
+    variable = T
+    other_variable = T_reconst
+  []
+[]
+
+[Debug]
+  show_var_residual_norms = true
+[]
+
+[Postprocessors]
+[]
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/2d.i b/modules/combined/examples/stochastic/laser_welding_dimred/2d.i
new file mode 100644
index 000000000000..cd6a8537c3af
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/2d.i
@@ -0,0 +1,52 @@
+!include parameters.i
+
+!include mesh.i
+
+[Variables]
+  [vel]
+    family = LAGRANGE_VEC
+  []
+  [T]
+  []
+  [p]
+  []
+  [disp_x]
+  []
+  [disp_y]
+  []
+[]
+
+!include physics_objects.i
+
+[Executioner]
+  type = Transient
+  end_time = ${endtime}
+  dtmin = 1e-10
+  dtmax = 1e-5
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu       NONZERO'
+  petsc_options = '-snes_converged_reason -ksp_converged_reason -options_left'
+  solve_type = 'NEWTON'
+  line_search = 'none'
+  nl_max_its = 5
+  l_max_its = 100
+  [TimeStepper]
+    type = IterationAdaptiveDT
+    optimal_iterations = 5
+    iteration_window = 1
+    dt = ${timestep}
+    linear_iteration_ratio = 1e6
+    growth_factor = 1.25
+  []
+[]
+
+[Reporters]
+  [solution_storage]
+    type = SolutionContainer
+    execute_on = 'FINAL'
+  []
+[]
+
+[Debug]
+  show_var_residual_norms = true
+[]
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/gold/error_to_plot.csv b/modules/combined/examples/stochastic/laser_welding_dimred/gold/error_to_plot.csv
new file mode 100644
index 000000000000..22c54d1e64b1
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/gold/error_to_plot.csv
@@ -0,0 +1,91 @@
+l2error
+0.0006183203894430075
+0.00031513281439315006
+0.0003046557003804341
+0.0005687926610983354
+0.0011695636499280098
+0.0024801093559574374
+0.0015842709716921742
+0.00047182449004882994
+0.00024218147315261865
+0.0005732994476687629
+0.00048185425767627294
+0.0005939617606608555
+0.00042919076573924976
+0.0007122785308927393
+0.0005516797745442211
+0.0005134402329385604
+0.00015381195177398318
+0.0017863174852499424
+0.00046959727881158887
+0.0007002480822981567
+0.0010516082914074783
+0.0020028017819068928
+0.00038504797137411544
+0.0016939107153106613
+0.00105861140658033
+0.00021662353233608698
+0.000715138695786721
+0.0011628690419034772
+0.0008397878209945399
+0.0009385317834876104
+0.00043060810567277627
+0.0013859987545918242
+0.0005116725931340964
+0.0006533789917551344
+0.000881415832914347
+0.0001856457685744164
+0.0013105208476763205
+0.00013023366709599208
+0.00151888284500697
+0.00051681013787522
+4.886151491773302e-
+0.00145967934675554
+0.00322130506505876
+0.00065269614540999
+0.00090565149002342
+0.011998030363012996
+0.00172482961399477
+0.00814709663479354
+0.00106699758206230
+0.00098731994057100
+0.00444282618703061
+0.00449629775447601
+0.00132091487401975
+0.00450876709096194
+0.00179141426485217
+0.010996130198897936
+0.00027840067583499
+0.00039411559492492
+0.00153263430337879
+0.005986536779903415
+0.001936200332505432
+0.00115491708844850
+0.009017710804415113
+0.00221946076557828
+0.007677715695900901
+0.00037850379053105
+0.00024378898844825
+0.00144422768915079
+0.00019160821391303
+0.009463951412202293
+0.001707714011530348
+0.00613848331400484
+0.015208431887830474
+0.002435591844646813
+0.004656097696630657
+0.015760722049044797
+0.009451338195770338
+0.00041603282055812
+0.00039079860784553
+0.00031512144006310
+0.00178612394676126
+0.00138977179586526
+0.00138831348975234
+0.00070426172150619
+0.00588496066404157
+0.01087780084733199
+0.00115516993602150
+0.00182441502097275
+0.00093793268808103
+0.0006539791112933185
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/gold/sv_to_print.csv b/modules/combined/examples/stochastic/laser_welding_dimred/gold/sv_to_print.csv
new file mode 100644
index 000000000000..3f40469d285c
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/gold/sv_to_print.csv
@@ -0,0 +1,9 @@
+i,sv
+1,540509.1951650645
+2,9737.187020905265
+3,2834.865343287165
+4,425.2316758454986
+5,185.21623877539616
+6,175.43364204675154
+7,156.0540415286562
+8,116.57686815340323
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/gold/test_json.json b/modules/combined/examples/stochastic/laser_welding_dimred/gold/test_json.json
new file mode 100644
index 000000000000..032454a10363
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/gold/test_json.json
@@ -0,0 +1,61 @@
+{
+    "reporters": {
+        "matrix": {
+            "type": "StochasticMatrix",
+            "values": {
+                "results:converged": {
+                    "row_begin": 0,
+                    "row_end": 4,
+                    "type": "std::vector<bool>"
+                },
+                "results:l2error:value": {
+                    "row_begin": 0,
+                    "row_end": 4,
+                    "type": "std::vector<double>"
+                },
+                "test_0": {
+                    "row_begin": 0,
+                    "row_end": 4,
+                    "type": "std::vector<double>"
+                },
+                "test_1": {
+                    "row_begin": 0,
+                    "row_end": 4,
+                    "type": "std::vector<double>"
+                }
+            }
+        }
+    },
+    "time_steps": [
+        {
+            "matrix": {
+                "results:converged": [
+                    true,
+                    true,
+                    true,
+                    true
+                ],
+                "results:l2error:value": [
+                    0.006140899524593516,
+                    0.010505678416121461,
+                    0.0015668821790721986,
+                    0.005425419862385741
+                ],
+                "test_0": [
+                    0.00013864103172310882,
+                    0.00012707154029110033,
+                    0.0001395556158341649,
+                    0.00013137278815868362
+                ],
+                "test_1": [
+                    68.32632033768421,
+                    63.19447251327895,
+                    62.98261633261899,
+                    61.66746839355559
+                ]
+            },
+            "time": 2.0,
+            "time_step": 2
+        }
+    ]
+}
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/mesh.i b/modules/combined/examples/stochastic/laser_welding_dimred/mesh.i
new file mode 100644
index 000000000000..7b44c4d8de90
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/mesh.i
@@ -0,0 +1,13 @@
+[Mesh]
+  [cmg]
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = ${xmin}
+    xmax = ${xmax}
+    ymin = ${fparse ymin}
+    ymax = 0
+    nx = 161
+    ny = 50
+  []
+  displacements = 'disp_x disp_y'
+[]
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/parameters.i b/modules/combined/examples/stochastic/laser_welding_dimred/parameters.i
new file mode 100644
index 000000000000..92abebdbd27d
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/parameters.i
@@ -0,0 +1,15 @@
+# Process parameters
+scanning_speed=1.0 # m/s
+power=74 # W (this is the effective power so multiplied by eta)
+R=1.25e-4 # m (this is the effective radius)
+
+# Geometric parameters
+thickness=1e-4 # m
+xmin=-0.1e-3 # m
+xmax=0.75e-3 # m
+ymin=${fparse -thickness}
+surfacetemp=300 # K (temperature at the other side of the plate)
+
+# Time stepping parameters
+endtime=4e-4 # s
+timestep=${fparse endtime/1000} # s
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/physics_objects.i b/modules/combined/examples/stochastic/laser_welding_dimred/physics_objects.i
new file mode 100644
index 000000000000..b50d312b7421
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/physics_objects.i
@@ -0,0 +1,206 @@
+[ICs]
+  [T]
+    type = FunctionIC
+    variable = T
+    function = '(${surfacetemp} - 300) / ${thickness} * y + ${surfacetemp}'
+  []
+[]
+
+[Kernels]
+  [disp_x]
+    type = Diffusion
+    variable = disp_x
+  []
+  [disp_y]
+    type = Diffusion
+    variable = disp_y
+  []
+  [mass]
+    type = INSADMass
+    variable = p
+    use_displaced_mesh = true
+  []
+  [mass_pspg]
+    type = INSADMassPSPG
+    variable = p
+    use_displaced_mesh = true
+  []
+  [momentum_time]
+    type = INSADMomentumTimeDerivative
+    variable = vel
+    use_displaced_mesh = true
+  []
+  [momentum_advection]
+    type = INSADMomentumAdvection
+    variable = vel
+    use_displaced_mesh = true
+  []
+  [momentum_mesh_advection]
+    type = INSADMomentumMeshAdvection
+    variable = vel
+    disp_x = disp_x
+    disp_y = disp_y
+    use_displaced_mesh = true
+  []
+  [momentum_viscous]
+    type = INSADMomentumViscous
+    variable = vel
+    use_displaced_mesh = true
+  []
+  [momentum_pressure]
+    type = INSADMomentumPressure
+    variable = vel
+    pressure = p
+    integrate_p_by_parts = true
+    use_displaced_mesh = true
+  []
+  [momentum_supg]
+    type = INSADMomentumSUPG
+    variable = vel
+    material_velocity = relative_velocity
+    use_displaced_mesh = true
+  []
+  [temperature_time]
+    type = INSADHeatConductionTimeDerivative
+    variable = T
+    use_displaced_mesh = true
+  []
+  [temperature_advection]
+    type = INSADEnergyAdvection
+    variable = T
+    use_displaced_mesh = true
+  []
+  [temperature_mesh_advection]
+    type = INSADEnergyMeshAdvection
+    variable = T
+    disp_x = disp_x
+    disp_y = disp_y
+    use_displaced_mesh = true
+  []
+  [temperature_conduction]
+    type = ADHeatConduction
+    variable = T
+    thermal_conductivity = 'k'
+    use_displaced_mesh = true
+  []
+  [temperature_supg]
+    type = INSADEnergySUPG
+    variable = T
+    velocity = vel
+    use_displaced_mesh = true
+  []
+[]
+
+[BCs]
+  [x_no_disp]
+    type = DirichletBC
+    variable = disp_x
+    boundary = 'bottom'
+    value = 0
+  []
+  [y_no_disp]
+    type = DirichletBC
+    variable = disp_y
+    boundary = 'bottom'
+    value = 0
+  []
+  [no_slip]
+    type = ADVectorFunctionDirichletBC
+    variable = vel
+    boundary = 'bottom right left'
+  []
+  [T_cold]
+    type = DirichletBC
+    variable = T
+    boundary = 'bottom'
+    value = 300
+  []
+  [radiation_flux]
+    type = FunctionRadiativeBC
+    variable = T
+    boundary = 'top'
+    emissivity_function = '1'
+    Tinfinity = 300
+    stefan_boltzmann_constant = 5.67e-8
+    use_displaced_mesh = true
+  []
+  [weld_flux]
+    type = GaussianEnergyFluxBC
+    variable = T
+    boundary = 'top'
+    P0 = ${power}
+    R = ${R}
+    x_beam_coord = '${scanning_speed}*t'
+    y_beam_coord = '0'
+    use_displaced_mesh = true
+  []
+  [vapor_recoil]
+    type = INSADVaporRecoilPressureMomentumFluxBC
+    variable = vel
+    boundary = 'top'
+    use_displaced_mesh = true
+  []
+  [surface_tension]
+    type = INSADSurfaceTensionBC
+    variable = vel
+    boundary = 'top'
+    use_displaced_mesh = true
+  []
+  [displace_x_top]
+    type = INSADDisplaceBoundaryBC
+    boundary = 'top'
+    variable = 'disp_x'
+    velocity = 'vel'
+    component = 0
+    associated_subdomain = 0
+  []
+  [displace_y_top]
+    type = INSADDisplaceBoundaryBC
+    boundary = 'top'
+    variable = 'disp_y'
+    velocity = 'vel'
+    component = 1
+    associated_subdomain = 0
+  []
+  [displace_x_top_dummy]
+    type = INSADDummyDisplaceBoundaryIntegratedBC
+    boundary = 'top'
+    variable = 'disp_x'
+    velocity = 'vel'
+    component = 0
+  []
+  [displace_y_top_dummy]
+    type = INSADDummyDisplaceBoundaryIntegratedBC
+    boundary = 'top'
+    variable = 'disp_y'
+    velocity = 'vel'
+    component = 1
+  []
+[]
+
+[Materials]
+  [ins_mat]
+    type = INSADStabilized3Eqn
+    velocity = vel
+    pressure = p
+    temperature = T
+    use_displaced_mesh = true
+  []
+  [steel]
+    type = LaserWeld316LStainlessSteel
+    temperature = T
+    use_displaced_mesh = true
+  []
+  [steel_boundary]
+    type = LaserWeld316LStainlessSteelBoundary
+    boundary = 'top'
+    temperature = T
+    use_displaced_mesh = true
+  []
+  [const]
+    type = GenericConstantMaterial
+    prop_names = 'abs sb_constant'
+    prop_values = '1 5.67e-8'
+    use_displaced_mesh = true
+  []
+[]
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/plot_sv.py b/modules/combined/examples/stochastic/laser_welding_dimred/plot_sv.py
new file mode 100644
index 000000000000..4b95f57210a6
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/plot_sv.py
@@ -0,0 +1,94 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import json
+plt.rcParams.update({
+    "text.usetex": True,
+    "font.family": "sans-serif",
+    "font.sans-serif": "Helvetica",
+})
+
+data = {}
+with open("pod_mapping_train_out.json", 'r') as file:
+  # Load the JSON data
+  data = json.load(file)
+
+aux_sv = data["time_steps"][0]["svd_output_aux"]["singular_values"]
+sol_sv = data["time_steps"][0]["svd_output_sol"]["singular_values"]
+
+colors = ['b','g','k']
+markers = ['^','o','v']
+aux_names = [r'$\mathrm{Velocity~x}$', r'$\mathrm{Velocity~y}$']
+aux_data = [aux_sv["vel_x_aux"],aux_sv["vel_y_aux"]]
+sum_aux_data = [sum([j*j for j in aux_data[i]]) for i in range(len(aux_data))]
+for i in range(len(sum_aux_data)):
+  aux_data[i] = [j*j/sum_aux_data[i] for j in aux_data[i]]
+
+sol_names = [r'$\mathrm{Temperature}$', r'$\mathrm{Displacement~x}$', r'$\mathrm{Displacement~y}$']
+sol_data = [sol_sv["T"],sol_sv["disp_x"],sol_sv["disp_y"]]
+sum_sol_data = [sum([j*j for j in sol_data[i]]) for i in range(len(sol_data))]
+sol_data = [i for i in sol_data]
+for i in range(len(sum_sol_data)):
+  sol_data[i] = [j*j/sum_sol_data[i] for j in sol_data[i]]
+
+# Create a figure and axis object
+fig, ax = plt.subplots(figsize=(8, 6))
+x = [i for i in range(1,21)]
+
+for i in range(len(aux_data)):
+  # Plotting the singular values
+  ax.semilogy(x,aux_data[i][0:20], marker=markers[i], linestyle='-', color=colors[i],
+          markersize=8, linewidth=2, label=aux_names[i])
+
+# Set custom axis ranges
+ax.set_xlim(1, 20)
+ax.set_ylim(1e-8, 1)
+ax.set_xticks(x)
+ax.tick_params(axis='y', which='both')
+ax.tick_params(axis='both', which='major', labelsize=20)
+
+ax.set_xlabel(r'$\mathrm{Index}$', fontsize=20)
+ax.set_ylabel(r'$\mathrm{Normalized~Squared~Singular~Value}$', fontsize=20)
+
+
+# Set grid
+ax.grid(True, linestyle='--', alpha=0.5)
+
+# Add legend
+ax.legend(fontsize=20, loc='best')
+
+# Tight layout
+fig.tight_layout()
+
+# Save plot to PDF
+plt.savefig('vel_aux.pdf')
+
+# Create a figure and axis object
+fig, ax = plt.subplots(figsize=(8, 6))
+
+for i in range(len(sol_data)):
+  # Plotting the singular values
+  ax.semilogy(x,sol_data[i][0:20], marker=markers[i], linestyle='-', color=colors[i],
+          markersize=8, linewidth=2, label=sol_names[i])
+
+# Set custom axis ranges
+ax.set_xlim(1, 20)
+ax.set_ylim(1e-10, 1)
+ax.set_xticks(x)
+ax.tick_params(axis='y', which='both')
+ax.tick_params(axis='both', which='major', labelsize=20)
+
+ax.set_xlabel(r'$\mathrm{Index}$', fontsize=20)
+ax.set_ylabel(r'$\mathrm{Normalized~Squared~Singular~Value}$', fontsize=20)
+
+
+# Set grid
+ax.grid(True, linestyle='--', alpha=0.5)
+
+# Add legend
+ax.legend(fontsize=20, loc='best')
+
+# Tight layout
+fig.tight_layout()
+
+# Save plot to PDF
+plt.savefig('sol_aux.pdf')
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/test.i b/modules/combined/examples/stochastic/laser_welding_dimred/test.i
new file mode 100644
index 000000000000..58d8d4a41d6e
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/test.i
@@ -0,0 +1,73 @@
+[StochasticTools]
+[]
+
+[Distributions]
+  [R_dist]
+    type = Uniform
+    lower_bound = 1.25E-4
+    upper_bound = 1.55E-4
+  []
+  [power_dist]
+    type = Uniform
+    lower_bound = 60
+    upper_bound = 70
+  []
+[]
+
+[Samplers]
+  [test]
+    type = MonteCarlo
+    num_rows = 90
+    distributions = 'R_dist power_dist'
+    execute_on = PRE_MULTIAPP_SETUP
+    min_procs_per_row = 2
+    max_procs_per_row = 2
+    seed=42
+  []
+[]
+
+[MultiApps]
+  [worker]
+    type = SamplerFullSolveMultiApp
+    input_files = 2d-reconst.i
+    sampler = test
+    mode = batch-reset
+    min_procs_per_app = 2
+    max_procs_per_app = 2
+  []
+[]
+
+[Controls]
+  [cmdline]
+    type = MultiAppSamplerControl
+    multi_app = worker
+    sampler = test
+    param_names = 'R power'
+  []
+[]
+
+[Transfers]
+  [results]
+    type = SamplerReporterTransfer
+    from_multi_app = worker
+    sampler = test
+    stochastic_reporter = matrix
+    from_reporter = 'l2error/value'
+  []
+[]
+
+[Reporters]
+  [matrix]
+    type = StochasticMatrix
+    sampler = test
+    parallel_type = ROOT
+  []
+[]
+
+[Outputs]
+  [json]
+    type = JSON
+    execute_on = FINAL
+    execute_system_information_on = NONE
+  []
+[]
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/tests b/modules/combined/examples/stochastic/laser_welding_dimred/tests
new file mode 100644
index 000000000000..7adb1f75735c
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/tests
@@ -0,0 +1,40 @@
+[Tests]
+  design = '/stochastic_tools/index.md'
+  issues = '#28200'
+
+  [train]
+    type = CheckFiles
+    input = train.i
+    cli_args = 'Samplers/train/num_rows=4
+                worker:/Mesh/cmg/nx=30
+                worker:/Mesh/cmg/ny=5
+                worker:endtime=1e-4
+                worker:timestep=1e-6
+                VariableMappings/pod/num_modes_to_compute=4
+                Covariance/lmc/num_outputs=4'
+    check_files = 'train_mogp_out_mogp.rd/data
+                   train_pod_out_pod.rd/data
+                   train_json.json'
+    requirement = "The system shall be able to train a surrogate model on the coefficients obtained by dimensionality reduction."
+    heavy = true
+    min_parallel = 8
+    max_parallel = 8
+    allow_test_objects = true
+  []
+  [test]
+    type = JSONDiff
+    input = test.i
+    cli_args = 'Samplers/test/num_rows=4
+                worker:/Mesh/cmg/nx=30
+                worker:/Mesh/cmg/ny=5
+                worker:endtime=1e-4
+                worker:timestep=1e-6'
+    jsondiff = 'test_json.json'
+    requirement = "The system shall be able to test the full field reconstruction using inverse mapping."
+    heavy = true
+    min_parallel = 8
+    max_parallel = 8
+    allow_test_objects = true
+    prereq = train
+  []
+[]
diff --git a/modules/combined/examples/stochastic/laser_welding_dimred/train.i b/modules/combined/examples/stochastic/laser_welding_dimred/train.i
new file mode 100644
index 000000000000..e1d454507224
--- /dev/null
+++ b/modules/combined/examples/stochastic/laser_welding_dimred/train.i
@@ -0,0 +1,150 @@
+[StochasticTools]
+[]
+
+[Distributions]
+  [R_dist]
+    type = Uniform
+    lower_bound = 1.25E-4
+    upper_bound = 1.55E-4
+  []
+  [power_dist]
+    type = Uniform
+    lower_bound = 60
+    upper_bound = 74
+  []
+[]
+
+[Samplers]
+  [train]
+    type = MonteCarlo
+    num_rows = 45
+    distributions = 'R_dist power_dist'
+    execute_on = PRE_MULTIAPP_SETUP
+    min_procs_per_row = 2
+    max_procs_per_row = 2
+  []
+[]
+
+[MultiApps]
+  [worker]
+    type = SamplerFullSolveMultiApp
+    input_files = 2d.i
+    sampler = train
+    mode = batch-reset
+    min_procs_per_app = 2
+    max_procs_per_app = 2
+  []
+[]
+
+[Controls]
+  [cmdline]
+    type = MultiAppSamplerControl
+    multi_app = worker
+    sampler = train
+    param_names = 'R power'
+  []
+[]
+
+[Transfers]
+  [solution_transfer]
+    type = SerializedSolutionTransfer
+    parallel_storage = parallel_storage
+    from_multi_app = worker
+    sampler = train
+    solution_container = solution_storage
+    variables = "T"
+    serialize_on_root = true
+  []
+[]
+
+[Reporters]
+  [parallel_storage]
+    type = ParallelSolutionStorage
+    variables = "T"
+    outputs = none
+  []
+  [reduced_sol]
+    type = MappingReporter
+    sampler = train
+    parallel_storage = parallel_storage
+    mapping = pod
+    variables = "T"
+    outputs = json
+    execute_on = final
+    parallel_type = ROOT
+  []
+  [matrix]
+    type = StochasticMatrix
+    sampler = train
+    parallel_type = ROOT
+  []
+  [svd]
+    type = SingularTripletReporter
+    pod_mapping = pod
+    variables = "T"
+    execute_on = final
+  []
+[]
+
+[VariableMappings]
+  [pod]
+    type = PODMapping
+    solution_storage = parallel_storage
+    variables = "T"
+    num_modes_to_compute = '8'
+    extra_slepc_options = "-svd_monitor_all"
+  []
+[]
+
+[Trainers]
+  [mogp]
+    type = GaussianProcessTrainer
+    execute_on = final
+    covariance_function = 'lmc'
+    standardize_params = 'true'
+    standardize_data = 'true'
+    sampler = train
+    response_type = vector_real
+    response = reduced_sol/T_pod
+    tune_parameters = 'lmc:acoeff_0 lmc:lambda_0 covar:signal_variance covar:length_factor'
+    tuning_min = '1e-9 1e-9 1e-9 1e-9'
+    tuning_max = '1e16 1e16 1e16  1e16'
+    num_iters = 10000
+    learning_rate = 0.0005
+    show_every_nth_iteration = 10
+  []
+[]
+
+[Covariance]
+  [covar]
+    type = SquaredExponentialCovariance
+    signal_variance = 1.0
+    noise_variance = 0.0
+    length_factor = '0.1 0.1'
+  []
+  [lmc]
+    type = LMC
+    covariance_functions = covar
+    num_outputs = 8
+    num_latent_funcs = 1
+  []
+[]
+
+
+[Outputs]
+  [json]
+    type = JSON
+    execute_on = FINAL
+    execute_system_information_on = NONE
+  []
+  [pod_out]
+    type = MappingOutput
+    mappings = pod
+    execute_on = FINAL
+  []
+  [mogp_out]
+    type = SurrogateTrainerOutput
+    trainers = mogp
+    execute_on = FINAL
+  []
+[]
diff --git a/modules/combined/examples/stochastic/gold/lhs_uniform_out.json b/modules/combined/examples/stochastic/thermomech/gold/lhs_uniform_out.json
similarity index 100%
rename from modules/combined/examples/stochastic/gold/lhs_uniform_out.json
rename to modules/combined/examples/stochastic/thermomech/gold/lhs_uniform_out.json
diff --git a/modules/combined/examples/stochastic/gold/poly_chaos_uniform_out.json b/modules/combined/examples/stochastic/thermomech/gold/poly_chaos_uniform_out.json
similarity index 100%
rename from modules/combined/examples/stochastic/gold/poly_chaos_uniform_out.json
rename to modules/combined/examples/stochastic/thermomech/gold/poly_chaos_uniform_out.json
diff --git a/modules/combined/examples/stochastic/graphite_ring_thermomechanics.i b/modules/combined/examples/stochastic/thermomech/graphite_ring_thermomechanics.i
similarity index 100%
rename from modules/combined/examples/stochastic/graphite_ring_thermomechanics.i
rename to modules/combined/examples/stochastic/thermomech/graphite_ring_thermomechanics.i
diff --git a/modules/combined/examples/stochastic/lhs_uniform.i b/modules/combined/examples/stochastic/thermomech/lhs_uniform.i
similarity index 100%
rename from modules/combined/examples/stochastic/lhs_uniform.i
rename to modules/combined/examples/stochastic/thermomech/lhs_uniform.i
diff --git a/modules/combined/examples/stochastic/poly_chaos_train_uniform.i b/modules/combined/examples/stochastic/thermomech/poly_chaos_train_uniform.i
similarity index 100%
rename from modules/combined/examples/stochastic/poly_chaos_train_uniform.i
rename to modules/combined/examples/stochastic/thermomech/poly_chaos_train_uniform.i
diff --git a/modules/combined/examples/stochastic/poly_chaos_uniform.i b/modules/combined/examples/stochastic/thermomech/poly_chaos_uniform.i
similarity index 100%
rename from modules/combined/examples/stochastic/poly_chaos_uniform.i
rename to modules/combined/examples/stochastic/thermomech/poly_chaos_uniform.i
diff --git a/modules/combined/examples/stochastic/tests b/modules/combined/examples/stochastic/thermomech/tests
similarity index 100%
rename from modules/combined/examples/stochastic/tests
rename to modules/combined/examples/stochastic/thermomech/tests
diff --git a/modules/navier_stokes/doc/content/bib/navier_stokes.bib b/modules/navier_stokes/doc/content/bib/navier_stokes.bib
index fa2302bb7a65..501cf753c71e 100644
--- a/modules/navier_stokes/doc/content/bib/navier_stokes.bib
+++ b/modules/navier_stokes/doc/content/bib/navier_stokes.bib
@@ -487,3 +487,42 @@ @article{menter1994two
   pages={1598--1605},
   year={1994}
 }
+
+@article{cairncross2000finite,
+  title={A finite element method for free surface flows of incompressible fluids in three dimensions. Part I. Boundary fitted mesh motion},
+  author={Cairncross, Richard A and Schunk, P Randall and Baer, Thomas A and Rao, Rekha R and Sackinger, Phillip A},
+  journal={International journal for numerical methods in fluids},
+  volume={33},
+  number={3},
+  pages={375--403},
+  year={2000},
+  publisher={Wiley Online Library}
+}
+
+@techreport{kim1975thermophysical,
+  title={Thermophysical properties of stainless steels},
+  author={Kim, Choong S},
+  year={1975},
+  institution={Argonne National Lab., Ill.(USA)}
+}
+
+@article{pichler2022surface,
+  title={Surface tension and thermal conductivity of NIST SRM 1155a (AISI 316L stainless steel)},
+  author={Pichler, Peter and Leitner, Thomas and Kaschnitz, Erhard and Rattenberger, Johannes and Pottlacher, Gernot},
+  journal={International Journal of Thermophysics},
+  volume={43},
+  number={5},
+  pages={66},
+  year={2022},
+  publisher={Springer}
+}
+
+@article{chen2021numerical,
+  title={Numerical analysis of double track formation for selective laser melting of 316L stainless steel},
+  author={Chen, Xuehui and Mu, Weihao and Xu, Xin and Liu, Wei and Huang, Lei and Li, Hao},
+  journal={Applied Physics A},
+  volume={127},
+  pages={1--13},
+  year={2021},
+  publisher={Springer}
+}
diff --git a/modules/navier_stokes/doc/content/modules/navier_stokes/laser_welding.md b/modules/navier_stokes/doc/content/modules/navier_stokes/laser_welding.md
index b972fa91294d..99d5a3a5ec44 100644
--- a/modules/navier_stokes/doc/content/modules/navier_stokes/laser_welding.md
+++ b/modules/navier_stokes/doc/content/modules/navier_stokes/laser_welding.md
@@ -32,13 +32,19 @@ flux, an outgoing radiative heat flux is modeled using [FunctionRadiativeBC.md]
 
 Evaporation of material from the liquefied material surface helps
 drive momentum changes at the surface of the condensed phase; this effect is incorporated via the
-[INSADVaporRecoilPressureMomentumFluxBC.md] object
+[INSADVaporRecoilPressureMomentumFluxBC.md] object:
 
 !listing modules/navier_stokes/examples/laser-welding/3d.i block=BCs/vapor_recoil
 
+The surface tension of the liquid also contributes to the
+forces that determine the deformation of the surface. This effect is
+added using the [INSADSurfaceTensionBC.md] object:
+
+!listing modules/navier_stokes/examples/laser-welding/3d.i block=BCs/surface_tension
+
 These surface momentum and velocity
 changes are then translated into mesh displacement
-through [INSADDisplaceBoundaryBC.md] objects
+through [INSADDisplaceBoundaryBC.md] objects.
 
 !listing modules/navier_stokes/examples/laser-welding/3d.i block=BCs/displace_x_top BCs/displace_y_top BCs/displace_z_top
 
diff --git a/modules/navier_stokes/doc/content/source/bcs/INSADSurfaceTensionBC.md b/modules/navier_stokes/doc/content/source/bcs/INSADSurfaceTensionBC.md
new file mode 100644
index 000000000000..7f95f2a2b5b9
--- /dev/null
+++ b/modules/navier_stokes/doc/content/source/bcs/INSADSurfaceTensionBC.md
@@ -0,0 +1,22 @@
+# INSADSurfaceTensionBC
+
+This boundary condition implements effects of surface tension on a free
+surface using the following expression:
+
+\begin{equation}
+\vec{n}\cdot \sigma_{liq} = -2 \mathcal{H} \sigma \vec{n} - \nabla_s \sigma \,,
+\end{equation}
+
+where $\mathcal{H}$ is the mean curvature of the surface, while
+$\sigma$ describes the suface tension. In this context, $\nabla_s = (I-\vec{n}\vec{n})\cdot \nabla$
+is the surface gradient operator. Parameter [!param](/BCs/INSADSurfaceTensionBC/include_gradient_terms)
+can be used to enable or disable the second term in the expression.
+Disabling the second term would disable the Marangoni effect and would decrease the
+surface deformations. This decreases the fidelity of the model, with an increased robustness.
+The model is based on the one discussed in [!cite](cairncross2000finite).
+
+!syntax parameters /BCs/INSADSurfaceTensionBC
+
+!syntax inputs /BCs/INSADSurfaceTensionBC
+
+!syntax children /BCs/INSADSurfaceTensionBC
diff --git a/modules/navier_stokes/examples/laser-welding/2d.i b/modules/navier_stokes/examples/laser-welding/2d.i
index 6a44a403371b..d7757119edbe 100644
--- a/modules/navier_stokes/examples/laser-welding/2d.i
+++ b/modules/navier_stokes/examples/laser-welding/2d.i
@@ -206,6 +206,13 @@ R=1.8257418583505537e-4 # m
     boundary = 'top'
     use_displaced_mesh = true
   []
+  [surface_tension]
+    type = INSADSurfaceTensionBC
+    variable = vel
+    boundary = 'top'
+    use_displaced_mesh = true
+    include_gradient_terms = true
+  []
   [displace_x_top]
     type = INSADDisplaceBoundaryBC
     boundary = 'top'
diff --git a/modules/navier_stokes/examples/laser-welding/3d.i b/modules/navier_stokes/examples/laser-welding/3d.i
index 70a140f2496a..e9dd9435368f 100644
--- a/modules/navier_stokes/examples/laser-welding/3d.i
+++ b/modules/navier_stokes/examples/laser-welding/3d.i
@@ -195,6 +195,12 @@ sb=5.67e-8
     boundary = 'front'
     use_displaced_mesh = true
   []
+  [surface_tension]
+    type = INSADSurfaceTensionBC
+    variable = vel
+    boundary = 'front'
+    use_displaced_mesh = true
+  []
   [displace_x_top]
     type = INSADDisplaceBoundaryBC
     boundary = 'front'
diff --git a/modules/navier_stokes/examples/laser-welding/gold/2d_curvature_out.e-s004 b/modules/navier_stokes/examples/laser-welding/gold/2d_curvature_out.e-s004
new file mode 100644
index 000000000000..af3a0908d9f1
Binary files /dev/null and b/modules/navier_stokes/examples/laser-welding/gold/2d_curvature_out.e-s004 differ
diff --git a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e
index 6f847be3e304..5d077d25e5d0 100644
Binary files a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e and b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e differ
diff --git a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s002 b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s002
index 5111e27a523f..ba3cfc1c0213 100644
Binary files a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s002 and b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s002 differ
diff --git a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s003 b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s003
index f571aa1fbcbb..a03d16825c11 100644
Binary files a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s003 and b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s003 differ
diff --git a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s004 b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s004
index f479065316fa..b0991cf7ee2a 100644
Binary files a/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s004 and b/modules/navier_stokes/examples/laser-welding/gold/2d_exodus.e-s004 differ
diff --git a/modules/navier_stokes/examples/laser-welding/tests b/modules/navier_stokes/examples/laser-welding/tests
index f125c6729d4d..401b545af954 100644
--- a/modules/navier_stokes/examples/laser-welding/tests
+++ b/modules/navier_stokes/examples/laser-welding/tests
@@ -1,6 +1,6 @@
 [Tests]
   design = 'laser_welding.md'
-  issues = '#24462'
+  issues = '#24462 #28200'
   [ALE]
     requirement = 'The system shall be able to simulate the rastering of a laser spot around an initially solid material which melts and then vaporizes'
     [2d-fe]
@@ -12,6 +12,21 @@
       valgrind = 'heavy'
       strumpack = true
       detail = 'in two dimensions, using a finite element discretization for all variables, and'
+      abs_zero = 1e-8
+      # PR #26848. Clang 16 Apple Si is not compatible.
+      machine = X86_64
+    []
+    [2d-fe-with-curvature]
+      input = 2d.i
+      exodiff = '2d_curvature_out.e-s004'
+      type = Exodiff
+      allow_test_objects = true
+      cli_args = "Executioner/num_steps=4 Mesh/elem_type=QUAD9 Variables/vel/order=SECOND Variables/disp_x/order=SECOND Variables/disp_y/order=SECOND  Postprocessors/active='' BCs/surface_tension/include_gradient_terms=false BCs/weld_flux/x_beam_coord=-0.35e-3 surfacetemp=2800 timestep=1e-6 Outputs/file_base=2d_curvature_out  Adaptivity/Markers/active='errorfrac_T combo' Adaptivity/Markers/combo/markers='errorfrac_T'"
+      valgrind = 'heavy'
+      strumpack = true
+      min_parallel = 8 # mesh refinement is very sensitive to number of cores here and before the 4th step we dont see deformation
+      max_parallel = 8 # mesh refinement is very sensitive to number of cores here and before the 4th step we dont see deformation
+      detail = 'in two dimensions, using a curvature-enabling second order mesh with finite element discretization, and'
       # PR #26848. Clang 16 Apple Si is not compatible.
       machine = X86_64
     []
diff --git a/modules/navier_stokes/include/bcs/INSADSurfaceTensionBC.h b/modules/navier_stokes/include/bcs/INSADSurfaceTensionBC.h
new file mode 100644
index 000000000000..8a3b3913a34f
--- /dev/null
+++ b/modules/navier_stokes/include/bcs/INSADSurfaceTensionBC.h
@@ -0,0 +1,44 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "ADVectorIntegratedBC.h"
+
+/**
+ * A class that imparts a surface tension on the momentum equation
+ * The treatment is based on:
+ * Cairncross RA, Schunk PR, Baer TA, Rao RR, Sackinger PA. A finite element method for free surface
+ * flows of incompressible fluids in three dimensions. Part I. Boundary fitted mesh motion.
+ * International journal for numerical methods in fluids. 2000 Jun 15;33(3):375-403.
+ */
+class INSADSurfaceTensionBC : public ADVectorIntegratedBC
+{
+public:
+  static InputParameters validParams();
+
+  INSADSurfaceTensionBC(const InputParameters & parameters);
+
+protected:
+  virtual ADReal computeQpResidual() override;
+
+  /// The surface tension terms
+  const ADMaterialProperty<RealVectorValue> & _surface_term_curvature;
+  const ADMaterialProperty<RealVectorValue> & _surface_term_gradient1;
+  const ADMaterialProperty<RealVectorValue> & _surface_term_gradient2;
+
+private:
+  /// If the surface tension should include the gradient terms
+  /// (increases fidelity, decreases stability)
+  const bool _include_gradient_terms;
+
+  /// Curvature force multiplier. The reason behind this is that
+  /// libmesh has a different sign convention for 2D and 3D.
+  const Real _curvature_factor;
+};
diff --git a/modules/navier_stokes/src/bcs/INSADSurfaceTensionBC.C b/modules/navier_stokes/src/bcs/INSADSurfaceTensionBC.C
new file mode 100644
index 000000000000..c8f58663a5a3
--- /dev/null
+++ b/modules/navier_stokes/src/bcs/INSADSurfaceTensionBC.C
@@ -0,0 +1,44 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "INSADSurfaceTensionBC.h"
+
+registerMooseObject("NavierStokesApp", INSADSurfaceTensionBC);
+
+InputParameters
+INSADSurfaceTensionBC::validParams()
+{
+  InputParameters params = ADVectorIntegratedBC::validParams();
+  params.addClassDescription("Surface tension stresses.");
+  params.addParam<bool>("include_gradient_terms",
+                        false,
+                        "If the surface tension should include the gradient terms (increases "
+                        "fidelity, decreases stability)");
+  return params;
+}
+
+INSADSurfaceTensionBC::INSADSurfaceTensionBC(const InputParameters & parameters)
+  : ADVectorIntegratedBC(parameters),
+    _surface_term_curvature(getADMaterialProperty<RealVectorValue>("surface_term_curvature")),
+    _surface_term_gradient1(getADMaterialProperty<RealVectorValue>("surface_term_gradient1")),
+    _surface_term_gradient2(getADMaterialProperty<RealVectorValue>("surface_term_gradient2")),
+    _include_gradient_terms(getParam<bool>("include_gradient_terms")),
+    _curvature_factor(_subproblem.mesh().dimension() == 3 ? 1.0 : -1.0)
+{
+}
+
+ADReal
+INSADSurfaceTensionBC::computeQpResidual()
+{
+  auto force = _curvature_factor * _surface_term_curvature[_qp];
+
+  if (_include_gradient_terms)
+    force += _surface_term_gradient1[_qp] + _surface_term_gradient2[_qp];
+  return -_test[_i][_qp] * force;
+}
diff --git a/modules/navier_stokes/test/include/materials/AriaLaserWeld304LStainlessSteelBoundary.h b/modules/navier_stokes/test/include/materials/AriaLaserWeld304LStainlessSteelBoundary.h
index a53b39af4301..904302c50c61 100644
--- a/modules/navier_stokes/test/include/materials/AriaLaserWeld304LStainlessSteelBoundary.h
+++ b/modules/navier_stokes/test/include/materials/AriaLaserWeld304LStainlessSteelBoundary.h
@@ -54,6 +54,8 @@ class AriaLaserWeld304LStainlessSteelBoundary : public Material
   ADMaterialProperty<Real> & _surface_tension;
   ADMaterialProperty<RealVectorValue> & _grad_surface_tension;
   const MooseArray<ADPoint> & _ad_normals;
+  const MooseArray<ADReal> & _ad_curvatures;
+  ADMaterialProperty<RealVectorValue> & _surface_term_curvature;
   ADMaterialProperty<RealVectorValue> & _surface_term_gradient1;
   ADMaterialProperty<RealVectorValue> & _surface_term_gradient2;
 };
diff --git a/modules/navier_stokes/test/include/materials/LaserWeld316LStainlessSteel.h b/modules/navier_stokes/test/include/materials/LaserWeld316LStainlessSteel.h
new file mode 100644
index 000000000000..68c51bfa818f
--- /dev/null
+++ b/modules/navier_stokes/test/include/materials/LaserWeld316LStainlessSteel.h
@@ -0,0 +1,99 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "Material.h"
+
+/**
+ * A material that computes 316L volumetric stainless steel properties relevant to doing laser
+ * welding modeling.
+ */
+class LaserWeld316LStainlessSteel : public Material
+{
+public:
+  static InputParameters validParams();
+
+  LaserWeld316LStainlessSteel(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpProperties();
+
+  /**
+   * An approximate for the dynamic viscosity is provided in paper:
+   * Kim, C.S., 1975. Thermophysical properties of stainless steels (No. ANL-75-55).
+   * Argonne National Lab., Ill.(USA)., (needed unit conversion)
+   */
+  const Real _c_mu0;
+  const Real _c_mu1;
+
+  /**
+   * The thermal conductivity is taken from:
+   * Pichler, Peter, et al. "Surface tension and thermal conductivity of NIST SRM 1155a
+   * (AISI 316L stainless steel)."
+   * International Journal of Thermophysics 43.5 (2022): 66.
+   */
+  const Real _c_k0_s;
+  const Real _c_k1_s;
+  const Real _c_k2_s;
+  const Real _c_k0_l;
+  const Real _c_k1_l;
+  const Real _c_k2_l;
+
+  /**
+   * An approximate for the speific heat is provided in paper:
+   * Kim, C.S., 1975. Thermophysical properties of stainless steels (No. ANL-75-55).
+   * Argonne National Lab., Ill.(USA)., (needed unit conversion)
+   */
+  const Real _c_cp0_s;
+  const Real _c_cp1_s;
+  const Real _c_cp0_l;
+
+  /**
+   * An approximate for the density is provided in paper:
+   * Kim, C.S., 1975. Thermophysical properties of stainless steels (No. ANL-75-55).
+   * Argonne National Lab., Ill.(USA)., (needed unit conversion)
+   */
+  const Real _c_rho0_s;
+  const Real _c_rho1_s;
+  const Real _c_rho2_s;
+  const Real _c_rho0_l;
+  const Real _c_rho1_l;
+  const Real _c_rho2_l;
+
+  /// Bounding temperature for the viscosity
+  const Real _Tmax;
+
+  /**
+   * The liquidus and solidus temperatures, both taken from:
+   * Pichler, Peter, et al. "Surface tension and thermal conductivity of NIST SRM 1155a
+   * (AISI 316L stainless steel)."
+   * International Journal of Thermophysics 43.5 (2022): 66.
+   */
+  const Real _Tl;
+  const Real _Ts;
+
+  /// We need to know the temperature and the gradient of the temperature
+  const ADVariableValue & _temperature;
+  const ADVariableGradient & _grad_temperature;
+
+  /// Material properties that get updated in this object
+  ADMaterialProperty<Real> & _mu;
+  ADMaterialProperty<Real> & _k;
+  ADMaterialProperty<Real> & _cp;
+  ADMaterialProperty<Real> & _rho;
+
+  /// The gradient of the thermal conductivity
+  ADMaterialProperty<RealVectorValue> & _grad_k;
+
+private:
+
+  /// If a constant density should be used (would discard buoyancy effects)
+  const bool _use_constant_density;
+};
diff --git a/modules/navier_stokes/test/include/materials/LaserWeld316LStainlessSteelBoundary.h b/modules/navier_stokes/test/include/materials/LaserWeld316LStainlessSteelBoundary.h
new file mode 100644
index 000000000000..7360565645cb
--- /dev/null
+++ b/modules/navier_stokes/test/include/materials/LaserWeld316LStainlessSteelBoundary.h
@@ -0,0 +1,64 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "Material.h"
+
+/**
+ * A material that computes 316L surface stainless steel boundary properties relevant to doing laser
+ * welding modeling (recoil pressure, surface tension).
+ */
+class LaserWeld316LStainlessSteelBoundary : public Material
+{
+public:
+  static InputParameters validParams();
+
+  LaserWeld316LStainlessSteelBoundary(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpProperties();
+
+  /**
+   * The vapor recoil pressure is taken from:
+   * Chen, Xuehui, et al. "Numerical analysis of double track formation for selective
+   * laser melting of 316L stainless steel." Applied Physics A 127 (2021): 1-13.
+   */
+  const Real _P0;
+  const Real _L_v;
+  const Real _M;
+  const Real _T_v;
+  const Real _R;
+
+  /**
+   * The surface tension is taken from:
+   * Pichler, Peter, et al. "Surface tension and thermal conductivity of NIST SRM 1155a
+   * (AISI 316L stainless steel)."
+   * International Journal of Thermophysics 43.5 (2022): 66.
+   */
+  const Real _c_gamma0;
+  const Real _c_gamma1;
+
+  /// The liquidus temperature
+  const Real _Tl;
+
+  /// Declaring the material properties
+  ADMaterialProperty<Real> & _rc_pressure;
+  ADMaterialProperty<Real> & _surface_tension;
+  ADMaterialProperty<RealVectorValue> & _grad_surface_tension;
+  ADMaterialProperty<RealVectorValue> & _surface_term_curvature;
+  ADMaterialProperty<RealVectorValue> & _surface_term_gradient1;
+  ADMaterialProperty<RealVectorValue> & _surface_term_gradient2;
+  const MooseArray<ADPoint> & _ad_normals;
+  const MooseArray<ADReal> & _ad_curvatures;
+
+  /// We need to know the temperature and the gradient of the temperature
+  const ADVariableValue & _temperature;
+  const ADVariableGradient & _grad_temperature;
+};
diff --git a/modules/navier_stokes/test/src/materials/AriaLaserWeld304LStainlessSteelBoundary.C b/modules/navier_stokes/test/src/materials/AriaLaserWeld304LStainlessSteelBoundary.C
index dc05524a5873..9beeb71ec106 100644
--- a/modules/navier_stokes/test/src/materials/AriaLaserWeld304LStainlessSteelBoundary.C
+++ b/modules/navier_stokes/test/src/materials/AriaLaserWeld304LStainlessSteelBoundary.C
@@ -66,6 +66,8 @@ AriaLaserWeld304LStainlessSteelBoundary::AriaLaserWeld304LStainlessSteelBoundary
     _grad_surface_tension(declareADProperty<RealVectorValue>(
         getParam<MaterialPropertyName>("grad_surface_tension_name"))),
     _ad_normals(_assembly.adNormals()),
+    _ad_curvatures(_assembly.adCurvatures()),
+    _surface_term_curvature(declareADProperty<RealVectorValue>("surface_term_curvature")),
     _surface_term_gradient1(declareADProperty<RealVectorValue>("surface_term_gradient1")),
     _surface_term_gradient2(declareADProperty<RealVectorValue>("surface_term_gradient2"))
 {
@@ -84,6 +86,8 @@ AriaLaserWeld304LStainlessSteelBoundary::computeQpProperties()
 
   _surface_tension[_qp] = _sigma0 + _alpha * (_temperature[_qp] - _T0);
   _grad_surface_tension[_qp] = _alpha * _grad_temperature[_qp];
+  _surface_term_curvature[_qp] =
+      -2. * _ad_curvatures[_qp] * _surface_tension[_qp] * _ad_normals[_qp];
   _surface_term_gradient1[_qp] = -_grad_surface_tension[_qp];
   _surface_term_gradient2[_qp] = _ad_normals[_qp] * (_ad_normals[_qp] * _grad_surface_tension[_qp]);
 }
diff --git a/modules/navier_stokes/test/src/materials/LaserWeld316LStainlessSteel.C b/modules/navier_stokes/test/src/materials/LaserWeld316LStainlessSteel.C
new file mode 100644
index 000000000000..60a8b84a49bc
--- /dev/null
+++ b/modules/navier_stokes/test/src/materials/LaserWeld316LStainlessSteel.C
@@ -0,0 +1,163 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "LaserWeld316LStainlessSteel.h"
+
+registerMooseObject("NavierStokesTestApp", LaserWeld316LStainlessSteel);
+
+InputParameters
+LaserWeld316LStainlessSteel::validParams()
+{
+  InputParameters params = Material::validParams();
+  params.addRequiredCoupledVar("temperature", "The temperature in K");
+
+  params.addParam<Real>("c_mu0", 2385.2, "mu0 coefficient in exp(m0/T-mu1)");
+  params.addParam<Real>("c_mu1", 0.5958, "mu1 coefficient in exp(m0/T-mu1)");
+
+  params.addParam<Real>("Tmax", 4000, "The maximum temperature for limiting viscosity");
+  params.addParam<Real>("Tl", 1708, "The liquidus temperature");
+  params.addParam<Real>("Ts", 1675, "The solidus temperature");
+
+  params.addParam<Real>("c_k0_s", 6.38, "k0 coefficient for the solid in k0+k1*T");
+  params.addParam<Real>("c_k1_s", 1.9E-2, "k1 coefficient for the solid in k0+k1*T");
+  params.addParam<Real>("c_k2_s", -2.45E-6, "k1 coefficient for the solid in k0+k1*T");
+  params.addParam<Real>("c_k0_l", 2.27, "k0 coefficient for the liquid in k0+k1*T");
+  params.addParam<Real>("c_k1_l", 1.76E-2, "k1 coefficien for the liquid in k0+k1*Tt");
+  params.addParam<Real>("c_k2_l", -1.39E-6, "k1 coefficient for the solid in k0+k1*T");
+
+  params.addParam<Real>("c_cp0_s", 459.29, "cp0 coefficient for the solid in cp0+cp1*T");
+  params.addParam<Real>("c_cp1_s", 0.1329, "cp1 coefficient for the solid in cp0+cp1*T");
+  params.addParam<Real>("c_cp0_l", 795.49, "cp0 coefficient for the liquid in cp0");
+
+  params.addParam<Real>(
+      "c_rho0_s", 8084.2, "The rho0 density for the solid in rho0+rho1*T+rho2*T*T");
+  params.addParam<Real>(
+      "c_rho1_s", -0.42086, "The rho1 density for the solid in rho0+rho1*T+rho2*T*T");
+  params.addParam<Real>(
+      "c_rho2_s", -3.8942E-5, "The rho2 density for the solid in rho0+rho1*T+rho2*T*T");
+  params.addParam<Real>(
+      "c_rho0_l", 7432.7, "The rho0 density for the liquid in rho0+rho1*T+rho2*T*T");
+  params.addParam<Real>(
+      "c_rho1_l", 0.039338, "The rho1 density for the liquid in rho0+rho1*T+rho2*T*T");
+  params.addParam<Real>(
+      "c_rho2_l", -1.8007E-4, "The rho2 density for the liquid in rho0+rho1*T+rho2*T*T");
+
+  params.addParam<MaterialPropertyName>(
+      "mu_name", "mu", "The name of the viscosity material property");
+  params.addParam<MaterialPropertyName>("k_name", "k", "The name of the thermal conductivity");
+  params.addParam<MaterialPropertyName>("cp_name", "cp", "The name of the specific heat capacity");
+  params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
+
+  params.addParam<bool>("use_constant_density", true, "If a constant density should be used (would discard buoyancy effects)");
+
+  return params;
+}
+
+LaserWeld316LStainlessSteel::LaserWeld316LStainlessSteel(
+    const InputParameters & parameters)
+  : Material(parameters),
+    _c_mu0(getParam<Real>("c_mu0")),
+    _c_mu1(getParam<Real>("c_mu1")),
+    _c_k0_s(getParam<Real>("c_k0_s")),
+    _c_k1_s(getParam<Real>("c_k1_s")),
+    _c_k2_s(getParam<Real>("c_k2_s")),
+    _c_k0_l(getParam<Real>("c_k0_l")),
+    _c_k1_l(getParam<Real>("c_k1_l")),
+    _c_k2_l(getParam<Real>("c_k2_l")),
+    _c_cp0_s(getParam<Real>("c_cp0_s")),
+    _c_cp1_s(getParam<Real>("c_cp1_s")),
+    _c_cp0_l(getParam<Real>("c_cp0_l")),
+    _c_rho0_s(getParam<Real>("c_rho0_s")),
+    _c_rho1_s(getParam<Real>("c_rho1_s")),
+    _c_rho2_s(getParam<Real>("c_rho2_s")),
+    _c_rho0_l(getParam<Real>("c_rho0_l")),
+    _c_rho1_l(getParam<Real>("c_rho1_l")),
+    _c_rho2_l(getParam<Real>("c_rho2_l")),
+    _Tmax(getParam<Real>("Tmax")),
+    _Tl(getParam<Real>("Tl")),
+    _Ts(getParam<Real>("Ts")),
+    _temperature(adCoupledValue("temperature")),
+    _grad_temperature(adCoupledGradient("temperature")),
+    _mu(declareADProperty<Real>(getParam<MaterialPropertyName>("mu_name"))),
+    _k(declareADProperty<Real>(getParam<MaterialPropertyName>("k_name"))),
+    _cp(declareADProperty<Real>(getParam<MaterialPropertyName>("cp_name"))),
+    _rho(declareADProperty<Real>(getParam<MaterialPropertyName>("rho_name"))),
+    _grad_k(declareADProperty<RealVectorValue>("grad_" + getParam<MaterialPropertyName>("k_name"))),
+    _use_constant_density(getParam<bool>("use_constant_density"))
+{
+}
+
+void
+LaserWeld316LStainlessSteel::computeQpProperties()
+{
+  ADReal That = _temperature[_qp] > _Tmax ? _Tmax : _temperature[_qp];
+  // First we check if it is enough to compute the solid properties only
+  if (_temperature[_qp] < _Ts)
+  {
+    // We are using an enormous viscosity for the solid phase to prevent it from moving.
+    // This approach was coined in:
+    // Noble, David R et al, Use of Aria to simulate laser weld pool dynamics for neutron generator production,
+    // 2007, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA
+    _mu[_qp] = 1e11 * 1e-3 * std::exp(_c_mu0 / That - _c_mu1);
+    _k[_qp] =
+        _c_k0_s + _c_k1_s * _temperature[_qp] + _c_k2_s * _temperature[_qp] * _temperature[_qp];
+    _grad_k[_qp] = _c_k1_s * _grad_temperature[_qp];
+    _cp[_qp] = _c_cp0_s + _c_cp1_s * _temperature[_qp];
+    _rho[_qp] = _c_rho0_s + _c_rho1_s * _temperature[_qp] +
+                _c_rho2_s * _temperature[_qp] * _temperature[_qp];
+  }
+  else
+  {
+    // This contains an artifically large lower viscosity at this point, to make sure
+    // the fluid simulations are stable. This viscosity is still below the one used in:
+    // Noble, David R et al, Use of Aria to simulate laser weld pool dynamics for neutron generator production,
+    // 2007, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA
+    const auto mu_l = std::max(0.030, 1e-3 * std::exp(_c_mu0 / That - _c_mu1));
+    const auto k_l =
+        _c_k0_l + _c_k1_l * _temperature[_qp] + _c_k2_l * _temperature[_qp] * _temperature[_qp];
+    const auto grad_k_l = _c_k1_l * _grad_temperature[_qp];
+    const auto cp_l = _c_cp0_l;
+    const auto rho_l = _c_rho0_l + _c_rho1_l * _temperature[_qp] +
+                       _c_rho2_l * _temperature[_qp] * _temperature[_qp];
+
+    // If we are between the solidus and liquidus temperatures we will need the
+    // solid properties too
+    if (_temperature[_qp] < _Tl && _temperature[_qp] > _Ts)
+    {
+      const auto liquid_fraction = (_temperature[_qp] - _Ts) / (_Tl - _Ts);
+
+      const auto mu_s = 1e11 * 1e-3 * std::exp(_c_mu0 / That - _c_mu1);
+      const auto k_s =
+          _c_k0_s + _c_k1_s * _temperature[_qp] + _c_k2_s * _temperature[_qp] * _temperature[_qp];
+      const auto grad_k_s = _c_k1_s * _grad_temperature[_qp];
+      const auto cp_s = _c_cp0_s + _c_cp1_s * _temperature[_qp];
+      const auto rho_s = _c_rho0_s + _c_rho1_s * _temperature[_qp] +
+                         _c_rho2_s * _temperature[_qp] * _temperature[_qp];
+
+      _mu[_qp] = liquid_fraction * mu_l + (1 - liquid_fraction) * mu_s;
+      _k[_qp] = liquid_fraction * k_l + (1 - liquid_fraction) * k_s;
+      _grad_k[_qp] = liquid_fraction * grad_k_l + (1 - liquid_fraction) * grad_k_s;
+      _cp[_qp] = liquid_fraction * cp_l + (1 - liquid_fraction) * cp_s;
+      _rho[_qp] = liquid_fraction * rho_l + (1 - liquid_fraction) * rho_s;
+    }
+    else
+    {
+      _mu[_qp] = mu_l;
+      _k[_qp] = k_l;
+      _grad_k[_qp] = grad_k_l;
+      _cp[_qp] = cp_l;
+      _rho[_qp] = rho_l;
+    }
+  }
+
+  // If we use constant density we override it with the density at the liquidus line
+  if (_use_constant_density)
+    _rho[_qp] = _c_rho0_l + _c_rho1_l * _Tl +
+                       _c_rho2_l * _Tl * _Tl;
+}
diff --git a/modules/navier_stokes/test/src/materials/LaserWeld316LStainlessSteelBoundary.C b/modules/navier_stokes/test/src/materials/LaserWeld316LStainlessSteelBoundary.C
new file mode 100644
index 000000000000..e6a04db26d5b
--- /dev/null
+++ b/modules/navier_stokes/test/src/materials/LaserWeld316LStainlessSteelBoundary.C
@@ -0,0 +1,92 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "LaserWeld316LStainlessSteelBoundary.h"
+#include "Assembly.h"
+
+registerMooseObject("NavierStokesTestApp", LaserWeld316LStainlessSteelBoundary);
+
+InputParameters
+LaserWeld316LStainlessSteelBoundary::validParams()
+{
+  InputParameters params = Material::validParams();
+  params.addParam<Real>("P0", 1e5, "The ambient pressure for the recoil pressure");
+  params.addParam<Real>("L_v", 7.45E6, "Latent heat from evaporation for the recoil pressure");
+  params.addParam<Real>("M", 56E-3, "Molar mass for the recoil pressure");
+  params.addParam<Real>("T_v", 3080, "The vaporization temperature for recoil pressure");
+  params.addParam<Real>("R", 8.314, "The gas constant for recoil pressure");
+  params.addParam<Real>("c_gamma0", 1.593, "Constant term in the surface tension");
+  params.addParam<Real>("c_gamma1", 1.9e-4, "Linear term multiplier in the surface tension");
+  params.addParam<Real>("Tl", 1708, "The liquidus temperature");
+  params.addRequiredCoupledVar("temperature", "The temperature in K");
+  params.addParam<MaterialPropertyName>("rc_pressure_name", "rc_pressure", "The recoil pressure");
+  params.addParam<MaterialPropertyName>(
+      "surface_tension_name", "surface_tension", "The surface tension");
+  params.addParam<MaterialPropertyName>(
+      "grad_surface_tension_name", "grad_surface_tension", "The gradient of the surface tension");
+  return params;
+}
+
+LaserWeld316LStainlessSteelBoundary::LaserWeld316LStainlessSteelBoundary(
+    const InputParameters & parameters)
+  : Material(parameters),
+    _P0(getParam<Real>("P0")),
+    _L_v(getParam<Real>("L_v")),
+    _M(getParam<Real>("M")),
+    _T_v(getParam<Real>("T_v")),
+    _R(getParam<Real>("R")),
+    _c_gamma0(getParam<Real>("c_gamma0")),
+    _c_gamma1(getParam<Real>("c_gamma1")),
+    _Tl(getParam<Real>("Tl")),
+    _rc_pressure(declareADProperty<Real>(getParam<MaterialPropertyName>("rc_pressure_name"))),
+    _surface_tension(
+        declareADProperty<Real>(getParam<MaterialPropertyName>("surface_tension_name"))),
+    _grad_surface_tension(declareADProperty<RealVectorValue>(
+        getParam<MaterialPropertyName>("grad_surface_tension_name"))),
+    _surface_term_curvature(declareADProperty<RealVectorValue>("surface_term_curvature")),
+    _surface_term_gradient1(declareADProperty<RealVectorValue>("surface_term_gradient1")),
+    _surface_term_gradient2(declareADProperty<RealVectorValue>("surface_term_gradient2")),
+    _ad_normals(_assembly.adNormals()),
+    _ad_curvatures(_assembly.adCurvatures()),
+    _temperature(adCoupledValue("temperature")),
+    _grad_temperature(adCoupledGradient("temperature"))
+{
+}
+
+void
+LaserWeld316LStainlessSteelBoundary::computeQpProperties()
+{
+  // Below the vaporization temperature we don't have recoil pressure
+  if (_temperature[_qp] < _T_v)
+    _rc_pressure[_qp] = 0;
+  else
+    _rc_pressure[_qp] =
+        0.54 * _P0 *
+        std::exp(_L_v * _M / _R * (_temperature[_qp] - _T_v) / (_temperature[_qp] * _T_v));
+
+  // The surface tension treatment methodology is from:
+  // Noble, David R et al, Use of Aria to simulate laser weld pool dynamics for neutron generator production,
+  // 2007, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA
+  //
+  // with the surface tension from:
+  //
+  // Pichler, Peter, et al. "Surface tension and thermal conductivity of NIST SRM 1155a
+  // (AISI 316L stainless steel)."
+  _surface_tension[_qp] = _c_gamma0 + _c_gamma1 * (_temperature[_qp] - _Tl);
+  _grad_surface_tension[_qp] = _c_gamma1 * _grad_temperature[_qp];
+
+  // These terms are needed for the treatment in:
+  // Cairncross RA, Schunk PR, Baer TA, Rao RR, Sackinger PA. A finite element method for free surface
+  // flows of incompressible fluids in three dimensions. Part I. Boundary fitted mesh motion.
+  // International journal for numerical methods in fluids. 2000 Jun 15;33(3):375-403.
+  _surface_term_curvature[_qp] =
+      -2. * _ad_curvatures[_qp] * _surface_tension[_qp] * _ad_normals[_qp];
+  _surface_term_gradient1[_qp] = -_grad_surface_tension[_qp];
+  _surface_term_gradient2[_qp] = _ad_normals[_qp] * (_ad_normals[_qp] * _grad_surface_tension[_qp]);
+}