diff --git a/simpegFLOW/Tests/test_Richards.py b/simpegFLOW/Tests/test_Richards.py
index 0e65ac5e..fd69adf4 100644
--- a/simpegFLOW/Tests/test_Richards.py
+++ b/simpegFLOW/Tests/test_Richards.py
@@ -3,6 +3,12 @@ from SimPEG import *
 from SimPEG.Tests.TestUtils import OrderTest, checkDerivative
 from scipy.sparse.linalg import dsolve
 from simpegFLOW import Richards
+try:
+    from pymatsolver import MumpsSolver
+    Solver = MumpsSolver
+except Exception, e:
+    pass
+
 
 TOL = 1E-8
 
@@ -76,6 +82,7 @@ class RichardsTests1D(unittest.TestCase):
         prob = Richards.RichardsProblem(M, mapping=E, timeSteps=[(40,3),(60,3)],
                                     boundaryConditions=bc, initialConditions=h,
                                     doNewton=False, method='mixed')
+        prob.Solver = Solver
 
         locs = np.r_[5.,10,15]
         times = prob.times[3:5]
@@ -147,6 +154,7 @@ class RichardsTests2D(unittest.TestCase):
         bc = np.r_[np.zeros(M.nCy*2),np.ones(M.nCx)*bc[0],np.ones(M.nCx)*bc[1]]
         h = np.zeros(M.nC) + bc[0]
         prob = Richards.RichardsProblem(M,E, timeSteps=[(40,3),(60,3)], boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed')
+        prob.Solver = Solver
 
         locs = Utils.ndgrid(np.array([5,7.]),np.array([5,15,25.]))
         times = prob.times[3:5]
@@ -202,5 +210,77 @@ class RichardsTests2D(unittest.TestCase):
         self.assertTrue(passed,True)
 
 
+
+class RichardsTests3D(unittest.TestCase):
+
+    def setUp(self):
+        M = Mesh.TensorMesh([np.ones(8),np.ones(20),np.ones(10)])
+
+        M.setCellGradBC(['neumann','neumann','dirichlet'])
+
+        params = Richards.Empirical.HaverkampParams().celia1990
+        params['Ks'] = np.log(params['Ks'])
+        E = Richards.Empirical.Haverkamp(M, **params)
+
+        bc = np.array([-61.5,-20.7])
+        bc = np.r_[np.zeros(M.nCy*M.nCz*2),np.zeros(M.nCx*M.nCz*2),np.ones(M.nCx*M.nCy)*bc[0],np.ones(M.nCx*M.nCy)*bc[1]]
+        h = np.zeros(M.nC) + bc[0]
+        prob = Richards.RichardsProblem(M,E, timeSteps=[(40,3),(60,3)], boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed')
+        prob.Solver = Solver
+
+        locs = Utils.ndgrid(np.r_[5,7.],np.r_[5,15.],np.r_[6,8.])
+        times = prob.times[3:5]
+        rxSat = Richards.RichardsRx(locs, times, 'saturation')
+        rxPre = Richards.RichardsRx(locs, times, 'pressureHead')
+        survey = Richards.RichardsSurvey([rxSat, rxPre])
+
+        prob.pair(survey)
+
+        self.h0 = h
+        self.M = M
+        self.Ks = params['Ks']
+        self.prob = prob
+        self.survey = survey
+
+    def test_Richards_getResidual_Newton(self):
+        self.prob.doNewton = True
+        m = self.Ks
+        passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False)
+        self.assertTrue(passed,True)
+
+    def test_Richards_getResidual_Picard(self):
+        self.prob.doNewton = False
+        m = self.Ks
+        passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False, expectedOrder=1)
+        self.assertTrue(passed,True)
+
+    def test_Adjoint(self):
+        v = np.random.rand(self.survey.nD)
+        z = np.random.rand(self.M.nC)
+        Hs = self.prob.fields(self.Ks)
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        tol = TOL*(10**int(np.log10(np.abs(zJv))))
+        passed = np.abs(vJz - zJv) < tol
+        print '3D: Richards Adjoint Test - PressureHead'
+        print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
+        self.assertTrue(passed,True)
+
+    def test_Sensitivity(self):
+        mTrue = self.Ks*np.ones(self.M.nC)
+        derChk = lambda m: [self.survey.dpred(m), lambda v: self.prob.Jvec(m, v)]
+        print '3D: Testing Richards Derivative'
+        passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
+        self.assertTrue(passed,True)
+
+    # def test_Sensitivity_full(self):
+    #     mTrue = self.Ks*np.ones(self.M.nC)
+    #     J = self.prob.Jfull(mTrue)
+    #     derChk = lambda m: [self.survey.dpred(m), J]
+    #     print '3D: Testing Richards Derivative FULL'
+    #     passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
+    #     self.assertTrue(passed,True)
+
+
 if __name__ == '__main__':
     unittest.main()