diff --git a/pinnicle/physics/__init__.py b/pinnicle/physics/__init__.py
index 95ef305..d46ace0 100644
--- a/pinnicle/physics/__init__.py
+++ b/pinnicle/physics/__init__.py
@@ -5,6 +5,7 @@
 from .continuity import *
 from .dummy import *
 from .stressbalance import *
+from .stressbalance_MC import *
 from .iceshelf import *
 from .timeinvariant import *
 from .friction import *
diff --git a/pinnicle/physics/constants.py b/pinnicle/physics/constants.py
index 69ec99f..e2a2f42 100644
--- a/pinnicle/physics/constants.py
+++ b/pinnicle/physics/constants.py
@@ -25,3 +25,9 @@ def __init__(self, **kwargs):
         self.variable_ub['u_base'] =  1.0e4/self.yts
         self.variable_lb['v_base'] = -1.0e4/self.yts
         self.variable_ub['v_base'] =  1.0e4/self.yts
+        self.variable_lb['D_dH'] = -1e3
+        self.variable_ub['D_dH'] = 1e3
+        self.variable_lb['D_smb'] = -1e3
+        self.variable_ub['D_smb'] = 1e3
+        self.variable_lb['R'] = -1e3
+        self.variable_ub['R'] = 1e3
diff --git a/pinnicle/physics/continuity.py b/pinnicle/physics/continuity.py
index f871d50..a38f46c 100644
--- a/pinnicle/physics/continuity.py
+++ b/pinnicle/physics/continuity.py
@@ -72,6 +72,80 @@ def _pde_jax(self, nn_input_var, nn_output_var): #{{{
         return self._pde(nn_input_var, nn_output_var) #}}}
     #}}}
 #}}}
+
+# Steady-state mass conservation (a=0) {{{
+class MCsteadyEquationParameter(EquationParameter, Constants):
+    """ default parameters for mass conservation
+    """
+    _EQUATION_TYPE = 'MC_steady' 
+    def __init__(self, param_dict={}):
+        # load necessary constants
+        Constants.__init__(self)
+        super().__init__(param_dict)
+
+    def set_default(self):
+        self.input = ['x', 'y']
+        self.output = ['u', 'v', 'H']
+        self.output_lb = [self.variable_lb[k] for k in self.output]
+        self.output_ub = [self.variable_ub[k] for k in self.output]
+        self.data_weights = [1.0e-8*self.yts**2, 1.0e-8*self.yts**2, 1.0e-6]
+        self.residuals = ["f"+self._EQUATION_TYPE]
+        self.pde_weights = [1.0e10]
+
+        # scalar variables: name:value
+        self.scalar_variables = {}
+class MC_steady(EquationBase): #{{{
+    """ MC on 2D problem with negligible smb and dH/dt
+    """
+    _EQUATION_TYPE = 'MC_steady' 
+    def __init__(self, parameters=MCEquationParameter()):
+        super().__init__(parameters)
+
+    def _pde(self, nn_input_var, nn_output_var): #{{{
+        """ residual of MC 2D PDE
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        # get the ids
+        xid = self.local_input_var["x"]
+        yid = self.local_input_var["y"]
+
+        uid = self.local_output_var["u"]
+        vid = self.local_output_var["v"]
+        Hid = self.local_output_var["H"]
+
+        # unpacking normalized output
+        u = slice_column(nn_output_var, uid)
+        v = slice_column(nn_output_var, vid)
+        H = slice_column(nn_output_var, Hid)
+
+        # spatial derivatives
+        u_x = jacobian(nn_output_var, nn_input_var, i=uid, j=xid)
+        H_x = jacobian(nn_output_var, nn_input_var, i=Hid, j=xid)
+        v_y = jacobian(nn_output_var, nn_input_var, i=vid, j=yid)
+        H_y = jacobian(nn_output_var, nn_input_var, i=Hid, j=yid)
+
+        # residual
+        f = H*u_x + H_x*u + H*v_y + H_y*v
+
+        return [f] #}}}
+    def _pde_jax(self, nn_input_var, nn_output_var): #{{{
+        """ residual of MC 2D PDE, jax version
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        return self._pde(nn_input_var, nn_output_var) #}}}
+    #}}}
+#}}}
+
+
+
+
+
 # time dependent mass conservation {{{
 class ThicknessEquationParameter(EquationParameter, Constants):
     """ default parameters for mass conservation
@@ -147,3 +221,121 @@ def _pde_jax(self, nn_input_var, nn_output_var): #{{{
         return self._pde(nn_input_var, nn_output_var) #}}}
     #}}}
 #}}}
+
+
+
+
+# D-HNN exact mass conservation {{{
+class MCexactEquationParameter(EquationParameter, Constants):
+    """ default parameters for mass conservation
+    """
+    _EQUATION_TYPE = 'MC_exact' 
+    def __init__(self, param_dict={}):
+        # load necessary constants
+        Constants.__init__(self)
+        super().__init__(param_dict)
+
+    def set_default(self):
+        self.input = ['x', 'y']
+        self.output = ['D_smb','D_dH','R', 'H']
+        self.output_lb = [self.variable_lb[k] for k in self.output]
+        self.output_ub = [self.variable_ub[k] for k in self.output]
+        self.data_weights = [1.0, 1.0, 1.0, 1.0e-6]
+        self.residuals = []
+        self.pde_weights = []
+
+        # scalar variables: name:value
+        self.scalar_variables = {}
+class MC_exact(EquationBase): #{{{
+    """ MC on 2D problem
+
+        for domains with negligible smb and dH/dt
+
+        u,v,a are defined based on two scalar fields D,R
+        in a way that automatically satisfies the MC
+    """
+    _EQUATION_TYPE = 'MC_exact' 
+    def __init__(self, parameters=MCexactEquationParameter()):
+        super().__init__(parameters)
+
+    def _pde(self, nn_input_var, nn_output_var): #{{{
+        """ no pde loss required
+            use data losses vel_mag_MC, u_MC, v_MC, a_MC
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        Hid = self.local_output_var["H"]
+        H = slice_column(nn_output_var, Hid)
+
+        return [] #}}}
+    
+    def _pde_jax(self, nn_input_var, nn_output_var): #{{{
+        """ residual of MC 2D PDE, jax version
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        return self._pde(nn_input_var, nn_output_var) #}}}
+    #}}}
+#}}}
+
+
+
+
+# D-HNN exact mass conservation {{{
+class MCSteadyexactEquationParameter(EquationParameter, Constants):
+    """ default parameters for mass conservation
+    """
+    _EQUATION_TYPE = 'MCSteady_exact' 
+    def __init__(self, param_dict={}):
+        # load necessary constants
+        Constants.__init__(self)
+        super().__init__(param_dict)
+
+    def set_default(self):
+        self.input = ['x', 'y']
+        self.output = ['R', 'H']
+        self.output_lb = [self.variable_lb[k] for k in self.output]
+        self.output_ub = [self.variable_ub[k] for k in self.output]
+        self.data_weights = [1.0, 1.0e-6]
+        self.residuals = []
+        self.pde_weights = []
+
+        # scalar variables: name:value
+        self.scalar_variables = {}
+class MCSteady_exact(EquationBase): #{{{
+    """ MC on 2D problem
+
+        for domains with negligible smb and dH/dt
+
+        u,v,a are defined based on two scalar fields D,R
+        in a way that automatically satisfies the MC
+    """
+    _EQUATION_TYPE = 'MCSteady_exact' 
+    def __init__(self, parameters=MCSteadyexactEquationParameter()):
+        super().__init__(parameters)
+
+    def _pde(self, nn_input_var, nn_output_var): #{{{
+        """ no pde loss required
+            use data losses vel_mag_MC, u_MC, v_MC, a_MC
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        return [] #}}}
+    
+    def _pde_jax(self, nn_input_var, nn_output_var): #{{{
+        """ residual of MC 2D PDE, jax version
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        return self._pde(nn_input_var, nn_output_var) #}}}
+    #}}}
+#}}}
+
diff --git a/pinnicle/physics/physics.py b/pinnicle/physics/physics.py
index fb157dc..575c4ea 100644
--- a/pinnicle/physics/physics.py
+++ b/pinnicle/physics/physics.py
@@ -125,3 +125,120 @@ def operator(self, pname):
             if p._EQUATION_TYPE.lower() == pname.lower():
                 return p.pde
 
+
+
+### transform functions for mass-conserving stressbalance:
+    def DR_xy(self, nn_input_var, nn_output_var):
+        """ transform D,R scalar fields of mass-conserving stressbalance
+            to u,v,a
+        """
+
+        xid = self.input_var.index('x')
+        yid = self.input_var.index('y')
+        
+        Rid = self.output_var.index('R')
+        R_x = jacobian(nn_output_var, nn_input_var, i=Rid, j=xid)
+        R_y = jacobian(nn_output_var, nn_input_var, i=Rid, j=yid)
+        
+        if "D_dH" in self.output_var:
+            # for MC_exact
+            Did = self.output_var.index('D_dH')
+            DdH_x = jacobian(nn_output_var, nn_input_var, i=Did, j=xid)
+            DdH_y = jacobian(nn_output_var, nn_input_var, i=Did, j=yid)
+        else:
+            # for MCSteady_exact
+            DdH_x = R_x*1e-32
+            DdH_y = R_y*1e-32
+
+        if "D_smb" in self.output_var:
+            # for MC_exact
+            Did = self.output_var.index('D_smb')
+            DdH_x = jacobian(nn_output_var, nn_input_var, i=Did, j=xid)
+            DdH_y = jacobian(nn_output_var, nn_input_var, i=Did, j=yid)
+        else:
+            # for MCSteady_exact
+            Dsmb_x = R_x*1e-32
+            Dsmb_y = R_y*1e-32
+
+        return [Dsmb_x, Dsmb_y, DdH_x, DdH_y, R_x, R_y]
+    
+    def DR_to_u(self, nn_input_var, nn_output_var):
+        """ recover u from scalar fields D,R
+        """
+        Hid = self.output_var.index('H')
+        H = slice_column(nn_output_var, Hid)
+        Dsmb_x, _, DdH_x, _, _, R_y = self.DR_xy(nn_input_var,nn_output_var)
+        # Hd = H.detach()
+        u = (Dsmb_x + DdH_x - R_y) / H 
+        return u
+    
+    def DR_to_v(self, nn_input_var, nn_output_var):
+        """ recover v from scalar fields D,R
+        """
+        Hid = self.output_var.index('H')
+        H = slice_column(nn_output_var, Hid)
+        _, Dsmb_y, _, DdH_y, R_x, _ = self.DR_xy(nn_input_var,nn_output_var)
+        # Hd = H.detach()
+        v = (Dsmb_y + DdH_y + R_x) / H # divide D_x by yts if dHdt in m/a
+        return v
+    
+    def DR_to_dH(self, nn_input_var, nn_output_var):
+        """ recover a from scalar fields D,R
+            a = div(grad(D))
+        """
+        xid = self.input_var.index('x')
+        yid = self.input_var.index('y')
+        _, _, DdH_x, DdH_y, _, _ = self.DR_xy(nn_input_var,nn_output_var)
+        DdH_xx = jacobian(DdH_x, nn_input_var, i=0, j=xid)
+        DdH_yy = jacobian(DdH_y, nn_input_var, i=0, j=yid)
+        dH = -1. * (DdH_xx + DdH_yy) ## == div(Hv)
+        return dH
+    
+    def DR_to_smb(self, nn_input_var, nn_output_var):
+        """ recover a from scalar fields D,R
+            a = div(grad(D))
+        """
+        xid = self.input_var.index('x')
+        yid = self.input_var.index('y')
+        Dsmb_x, Dsmb_y, _, _, _, _ = self.DR_xy(nn_input_var,nn_output_var)
+        Dsmb_xx = jacobian(Dsmb_x, nn_input_var, i=0, j=xid)
+        Dsmb_yy = jacobian(Dsmb_y, nn_input_var, i=0, j=yid)
+        smb = Dsmb_xx + Dsmb_yy ## == div(Hv)
+        return smb
+
+    def vel_mag_MC(self, nn_input_var, nn_output_var, X):
+        """ a wrapper for PointSetOperatorBC func call, Args need to follow the requirment by deepxde
+
+        Args: 
+            nn_input_var:  input tensor to the nn
+            nn_output_var: output tensor from the nn
+            X:  NumPy array of the collocation points defined on the boundary, required by deepxde
+        """
+        u = self.DR_to_u(nn_input_var,nn_output_var)
+        v = self.DR_to_v(nn_input_var,nn_output_var)
+        vel = ppow((bkd.square(u) + bkd.square(v) + 1.0e-30), 0.5)
+        return vel
+    
+    def dH_MC(self, nn_input_var, nn_output_var, X):
+        """ a wrapper for PointSetOperatorBC func call, Args need to follow the requirment by deepxde
+        """
+        dH = self.DR_to_dH(nn_input_var,nn_output_var)
+        return dH
+    
+    def smb_MC(self, nn_input_var, nn_output_var, X):
+        """ a wrapper for PointSetOperatorBC func call, Args need to follow the requirment by deepxde
+        """
+        smb = self.DR_to_smb(nn_input_var,nn_output_var)
+        return smb
+    
+    def u_MC(self, nn_input_var, nn_output_var, X):
+        """ a wrapper for PointSetOperatorBC func call, Args need to follow the requirment by deepxde
+        """
+        u = self.DR_to_u(nn_input_var,nn_output_var)
+        return u
+    
+    def v_MC(self, nn_input_var, nn_output_var, X):
+        """ a wrapper for PointSetOperatorBC func call, Args need to follow the requirment by deepxde
+        """
+        v = self.DR_to_v(nn_input_var,nn_output_var)
+        return v
\ No newline at end of file
diff --git a/pinnicle/physics/stressbalance_MC.py b/pinnicle/physics/stressbalance_MC.py
new file mode 100644
index 0000000..4ad175a
--- /dev/null
+++ b/pinnicle/physics/stressbalance_MC.py
@@ -0,0 +1,324 @@
+import deepxde as dde
+from deepxde.backend import jax, abs
+from . import EquationBase, Constants, Physics
+from ..parameter import EquationParameter
+from ..utils import slice_column, jacobian, slice_function_jax
+
+# mass-conserving SSA with taub
+class SSAMCTauEquationParameter(EquationParameter, Constants):
+    """ default parameters for SSA Taub
+    """
+    _EQUATION_TYPE = 'SSA_MC Taub' 
+    def __init__(self, param_dict={}):
+        # load necessary constants
+        Constants.__init__(self)
+        super().__init__(param_dict)
+
+    def set_default(self):
+        self.input = ['x', 'y']
+        self.output = ['D', 'R', 's', 'H', 'taub']
+        self.output_lb = [self.variable_lb[k] for k in self.output]
+        self.output_ub = [self.variable_ub[k] for k in self.output]
+        self.data_weights = [1.0, 1.0, 1.0e-6, 1.0e-6, 1.0e-10]#, 1.0e-2*self.yts**2]
+        self.residuals = ["f"+self._EQUATION_TYPE+"1", "f"+self._EQUATION_TYPE+"2"]
+        self.pde_weights = [1.0e-10, 1.0e-10]
+
+        # scalar variables: name:value
+        self.scalar_variables = {
+                'n': 3.0,               # exponent of Glen's flow law
+                'B':1.26802073401e+08   # -8 degree C, cuffey
+                }
+class SSA_MC_Taub(EquationBase): #{{{
+    """ SSA on 2D problem with uniform B, no friction law, but use taub=-beta*u
+    """
+    _EQUATION_TYPE = 'SSA_MC Taub' 
+    def __init__(self, parameters=SSAMCTauEquationParameter()):
+        super().__init__(parameters)
+    def _pde(self, nn_input_var, nn_output_var): #{{{
+        """ residual of SSA 2D PDEs
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        # get the ids
+        xid = self.local_input_var["x"]
+        yid = self.local_input_var["y"]
+
+        Did = self.local_output_var["D"]
+        Rid = self.local_output_var["R"]
+        sid = self.local_output_var["s"]
+        Hid = self.local_output_var["H"]
+        taubid = self.local_output_var["taub"]
+
+        # unpacking normalized output
+        H = slice_column(nn_output_var, Hid)
+        taub = slice_column(nn_output_var, taubid)
+        
+        # recovering u,v,a
+        D_x = jacobian(nn_output_var, nn_input_var, i=Did, j=xid)
+        D_y = jacobian(nn_output_var, nn_input_var, i=Did, j=yid)
+        R_x = jacobian(nn_output_var, nn_input_var, i=Rid, j=xid)
+        R_y = jacobian(nn_output_var, nn_input_var, i=Rid, j=yid)
+
+        # a = D_x + D_y ## == div(Hv)
+        u = (D_x - R_y) / H
+        v = (D_y + R_x) / H
+
+        # u = Physics.DR_to_u(self, nn_input_var, nn_output_var)
+        # v = Physics.DR_to_v(self, nn_input_var, nn_output_var)
+
+        # spatial derivatives
+        u_x = jacobian(u, nn_input_var, i=0, j=xid)
+        v_x = jacobian(v, nn_input_var, i=0, j=xid)
+        u_y = jacobian(u, nn_input_var, i=0, j=yid)
+        v_y = jacobian(v, nn_input_var, i=0, j=yid)
+        s_x = jacobian(nn_output_var, nn_input_var, i=sid, j=xid)
+        s_y = jacobian(nn_output_var, nn_input_var, i=sid, j=yid)
+
+        eta = 0.5*self.B *(u_x**2.0 + v_y**2.0 + 0.25*(u_y+v_x)**2.0 + u_x*v_y+self.eps)**(0.5*(1.0-self.n)/self.n)
+        # stress tensor
+        etaH = eta * H
+        B11 = etaH*(4*u_x + 2*v_y)
+        B22 = etaH*(4*v_y + 2*u_x)
+        B12 = etaH*(  u_y +   v_x)
+
+        # Getting the other derivatives
+        sigma11 = jacobian(B11, nn_input_var, i=0, j=xid)
+        sigma12 = jacobian(B12, nn_input_var, i=0, j=yid)
+
+        sigma21 = jacobian(B12, nn_input_var, i=0, j=xid)
+        sigma22 = jacobian(B22, nn_input_var, i=0, j=yid)
+
+
+        # compute the basal stress
+        u_norm = (u**2+v**2+self.eps**2)**0.5
+
+        f1 = sigma11 + sigma12 - abs(taub)*u/(u_norm) - self.rhoi*self.g*H*s_x
+        f2 = sigma21 + sigma22 - abs(taub)*v/(u_norm) - self.rhoi*self.g*H*s_y
+
+        return [f1, f2] #}}}
+    def _pde_jax(self, nn_input_var, nn_output_var): #{{{
+        """ residual of SSA 2D PDEs
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        # get the ids
+        xid = self.local_input_var["x"]
+        yid = self.local_input_var["y"]
+
+        Did = self.local_output_var["D"]
+        Rid = self.local_output_var["R"]
+        sid = self.local_output_var["s"]
+        Hid = self.local_output_var["H"]
+        taubid = self.local_output_var["taub"]
+
+        # unpacking normalized output
+        H = slice_column(nn_output_var, Hid)
+        taub = slice_column(nn_output_var, taubid)
+        
+        # recovering u,v,a
+        D_x = jacobian(nn_output_var, nn_input_var, i=Did, j=xid)
+        D_y = jacobian(nn_output_var, nn_input_var, i=Did, j=yid)
+        R_x = jacobian(nn_output_var, nn_input_var, i=Rid, j=xid)
+        R_y = jacobian(nn_output_var, nn_input_var, i=Rid, j=yid)
+
+        # a = D_x + D_y ## == div(Hv)
+        u = (D_x - R_y) / H
+        v = (D_y + R_x) / H
+
+        # get the spatial derivatives functions
+        u_x = jacobian(u, nn_input_var, i=0, j=xid, val=1)
+        v_x = jacobian(v, nn_input_var, i=0, j=xid, val=1)
+        u_y = jacobian(u, nn_input_var, i=0, j=yid, val=1)
+        v_y = jacobian(v, nn_input_var, i=0, j=yid, val=1)
+
+        # get variable function
+        H_func = lambda x: slice_function_jax(nn_output_var, x, Hid)
+        # stress tensor
+        etaH = lambda x: 0.5*H_func(x)*self.B *(u_x(x)**2.0 + v_y(x)**2.0 + 0.25*(u_y(x)+v_x(x))**2.0 + u_x(x)*v_y(x)+self.eps)**(0.5*(1.0-self.n)/self.n)
+
+        B11 = lambda x: etaH(x)*(4*u_x(x) + 2*v_y(x))
+        B22 = lambda x: etaH(x)*(4*v_y(x) + 2*u_x(x))
+        B12 = lambda x: etaH(x)*(  u_y(x) +   v_x(x))
+
+        # Getting the other derivatives
+        sigma11 = jacobian((jax.vmap(B11)(nn_input_var), B11), nn_input_var, i=0, j=xid)
+        sigma12 = jacobian((jax.vmap(B12)(nn_input_var), B12), nn_input_var, i=0, j=yid)
+
+        sigma21 = jacobian((jax.vmap(B12)(nn_input_var), B12), nn_input_var, i=0, j=xid)
+        sigma22 = jacobian((jax.vmap(B22)(nn_input_var), B22), nn_input_var, i=0, j=yid)
+
+        # compute the basal stress
+        s_x = jacobian(nn_output_var, nn_input_var, i=sid, j=xid)
+        s_y = jacobian(nn_output_var, nn_input_var, i=sid, j=yid)
+
+        u_norm = (u**2+v**2+self.eps**2)**0.5
+
+        f1 = sigma11 + sigma12 - abs(taub)*u/(u_norm) - self.rhoi*self.g*H*s_x
+        f2 = sigma21 + sigma22 - abs(taub)*v/(u_norm) - self.rhoi*self.g*H*s_y
+
+        return [f1, f2] #}}}
+    #}}}
+#}}}
+
+
+
+
+
+
+
+
+# mass-conserving SSA with taub
+# for regions with negligible smb and dH/dt
+class SSAMCsteadyTauEquationParameter(EquationParameter, Constants):
+    """ default parameters for SSA Taub
+    """
+    _EQUATION_TYPE = 'SSA_MCsteady Taub' 
+    def __init__(self, param_dict={}):
+        # load necessary constants
+        Constants.__init__(self)
+        super().__init__(param_dict)
+
+    def set_default(self):
+        self.input = ['x', 'y']
+        self.output = ['R', 's', 'H', 'taub']
+        self.output_lb = [self.variable_lb[k] for k in self.output]
+        self.output_ub = [self.variable_ub[k] for k in self.output]
+        self.data_weights = [1.0, 1.0e-5, 1.0e-5, 1.0e-10]
+        self.residuals = ["f"+self._EQUATION_TYPE+"1", "f"+self._EQUATION_TYPE+"2"]
+        self.pde_weights = [1.0e-10, 1.0e-10]
+
+        # scalar variables: name:value
+        self.scalar_variables = {
+                'n': 3.0,               # exponent of Glen's flow law
+                'B':1.26802073401e+08   # -8 degree C, cuffey
+                }
+class SSA_MCsteady_Taub(EquationBase): #{{{
+    """ SSA on 2D problem with uniform B, no friction law, but use taub
+    """
+    _EQUATION_TYPE = 'SSA_MCsteady Taub' 
+    def __init__(self, parameters=SSAMCsteadyTauEquationParameter()):
+        super().__init__(parameters)
+    def _pde(self, nn_input_var, nn_output_var): #{{{
+        """ residual of SSA 2D PDEs
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        # get the ids
+        xid = self.local_input_var["x"]
+        yid = self.local_input_var["y"]
+
+        Rid = self.local_output_var["R"]
+        sid = self.local_output_var["s"]
+        Hid = self.local_output_var["H"]
+        taubid = self.local_output_var["taub"]
+
+        # unpacking normalized output
+        H = slice_column(nn_output_var, Hid)
+        taub = slice_column(nn_output_var, taubid)
+        
+        # recovering u,v,a
+        R_x = jacobian(nn_output_var, nn_input_var, i=Rid, j=xid)
+        R_y = jacobian(nn_output_var, nn_input_var, i=Rid, j=yid)
+
+        # a = D_x + D_y ## == div(Hv)
+        u = -1. * R_y / H
+        v = R_x / H
+
+        # spatial derivatives
+        u_x = jacobian(u, nn_input_var, i=0, j=xid)
+        v_x = jacobian(v, nn_input_var, i=0, j=xid)
+        u_y = jacobian(u, nn_input_var, i=0, j=yid)
+        v_y = jacobian(v, nn_input_var, i=0, j=yid)
+        s_x = jacobian(nn_output_var, nn_input_var, i=sid, j=xid)
+        s_y = jacobian(nn_output_var, nn_input_var, i=sid, j=yid)
+
+        eta = 0.5*self.B *(u_x**2.0 + v_y**2.0 + 0.25*(u_y+v_x)**2.0 + u_x*v_y+self.eps)**(0.5*(1.0-self.n)/self.n)
+        # stress tensor
+        etaH = eta * H
+        B11 = etaH*(4*u_x + 2*v_y)
+        B22 = etaH*(4*v_y + 2*u_x)
+        B12 = etaH*(  u_y +   v_x)
+
+        # Getting the other derivatives
+        sigma11 = jacobian(B11, nn_input_var, i=0, j=xid)
+        sigma12 = jacobian(B12, nn_input_var, i=0, j=yid)
+
+        sigma21 = jacobian(B12, nn_input_var, i=0, j=xid)
+        sigma22 = jacobian(B22, nn_input_var, i=0, j=yid)
+
+
+        # compute the basal stress
+        u_norm = (u**2+v**2+self.eps**2)**0.5
+
+        f1 = sigma11 + sigma12 - abs(taub)*u/(u_norm) - self.rhoi*self.g*H*s_x
+        f2 = sigma21 + sigma22 - abs(taub)*v/(u_norm) - self.rhoi*self.g*H*s_y
+
+        return [f1, f2] #}}}
+    def _pde_jax(self, nn_input_var, nn_output_var): #{{{
+        """ residual of SSA 2D PDEs
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        # get the ids
+        xid = self.local_input_var["x"]
+        yid = self.local_input_var["y"]
+
+        Rid = self.local_output_var["R"]
+        sid = self.local_output_var["s"]
+        Hid = self.local_output_var["H"]
+        taubid = self.local_output_var["taub"]
+
+        # unpacking normalized output
+        H = slice_column(nn_output_var, Hid)
+        taub = slice_column(nn_output_var, taubid)
+        
+        # recovering u,v,a
+        R_x = jacobian(nn_output_var, nn_input_var, i=Rid, j=xid)
+        R_y = jacobian(nn_output_var, nn_input_var, i=Rid, j=yid)
+
+        # a = D_x + D_y ## == div(Hv)
+        u = -1. * R_y / H
+        v = R_x / H
+
+        # get the spatial derivatives functions
+        u_x = jacobian(u, nn_input_var, i=0, j=xid, val=1)
+        v_x = jacobian(v, nn_input_var, i=0, j=xid, val=1)
+        u_y = jacobian(u, nn_input_var, i=0, j=yid, val=1)
+        v_y = jacobian(v, nn_input_var, i=0, j=yid, val=1)
+
+        # get variable function
+        H_func = lambda x: slice_function_jax(nn_output_var, x, Hid)
+        # stress tensor
+        etaH = lambda x: 0.5*H_func(x)*self.B *(u_x(x)**2.0 + v_y(x)**2.0 + 0.25*(u_y(x)+v_x(x))**2.0 + u_x(x)*v_y(x)+self.eps)**(0.5*(1.0-self.n)/self.n)
+
+        B11 = lambda x: etaH(x)*(4*u_x(x) + 2*v_y(x))
+        B22 = lambda x: etaH(x)*(4*v_y(x) + 2*u_x(x))
+        B12 = lambda x: etaH(x)*(  u_y(x) +   v_x(x))
+
+        # Getting the other derivatives
+        sigma11 = jacobian((jax.vmap(B11)(nn_input_var), B11), nn_input_var, i=0, j=xid)
+        sigma12 = jacobian((jax.vmap(B12)(nn_input_var), B12), nn_input_var, i=0, j=yid)
+
+        sigma21 = jacobian((jax.vmap(B12)(nn_input_var), B12), nn_input_var, i=0, j=xid)
+        sigma22 = jacobian((jax.vmap(B22)(nn_input_var), B22), nn_input_var, i=0, j=yid)
+
+        # compute the basal stress
+        s_x = jacobian(nn_output_var, nn_input_var, i=sid, j=xid)
+        s_y = jacobian(nn_output_var, nn_input_var, i=sid, j=yid)
+
+        u_norm = (u**2+v**2+self.eps**2)**0.5
+
+        f1 = sigma11 + sigma12 - abs(taub)*u/(u_norm) - self.rhoi*self.g*H*s_x
+        f2 = sigma21 + sigma22 - abs(taub)*v/(u_norm) - self.rhoi*self.g*H*s_y
+
+        return [f1, f2] #}}}
+    #}}}
+#}}}
\ No newline at end of file
diff --git a/pinnicle/pinn.py b/pinnicle/pinn.py
index 9efc6d3..0229f77 100644
--- a/pinnicle/pinn.py
+++ b/pinnicle/pinn.py
@@ -266,6 +266,21 @@ def min_or_none(a, b):
                     elif d == "vel":
                         training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.vel_mag,
                                                                          batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
+                    elif d == "vel_MC":
+                        training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.vel_mag_MC,
+                                                                         batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
+                    elif d == "dH_MC":
+                        training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.dH_MC,
+                                                                         batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
+                    elif d == "smb_MC":
+                        training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.smb_MC,
+                                                                         batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
+                    elif d == "u_MC":
+                        training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.u_MC,
+                                                                         batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
+                    elif d == "v_MC":
+                        training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.v_MC,
+                                                                         batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
                     elif d == "sx":
                         training_temp.append(dde.icbc.PointSetOperatorBC(training_data.X[d], training_data.sol[d], self.physics.user_defined_gradient('s','x'),
                                                                          batch_size=min_or_none(self.params.training.mini_batch, training_data.sol[d].shape[0]), shuffle=True))
diff --git a/tests/test_pde.py b/tests/test_pde.py
index 0617078..235fa31 100644
--- a/tests/test_pde.py
+++ b/tests/test_pde.py
@@ -89,6 +89,16 @@ def test_SSA_Taub_pde_function():
     assert y[0].shape == (10,1)
     assert y[1].shape == (10,1)
 
+def test_SSA_MC_Taub_pde_function():
+    hp_local = dict(hp)
+    hp_local["equations"] = {"SSA_MC Taub":{}}
+    experiment = pinn.PINN(params=hp_local)
+    experiment.compile()
+    y = experiment.model.predict(experiment.model_data.X['u'], operator=experiment.physics.operator("SSA_MC Taub"))
+    assert len(y) == 2
+    assert y[0].shape == (10,1)
+    assert y[1].shape == (10,1)
+
 @pytest.mark.skipif(backend_name=="jax", reason="MOLHO is not implemented for jax")
 def test_MOLHO_pde_function():
     hp_local = dict(hp)