diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py
new file mode 100644
index 00000000..a86f9697
--- /dev/null
+++ b/SimPEG/forward/Richards.py
@@ -0,0 +1,219 @@
+from SimPEG.forward import Problem
+import numpy as np
+from SimPEG.utils import sdiag, mkvc, setKwargs
+
+
+class RichardsProblem(Problem):
+    """docstring for RichardsProblem"""
+
+    dt = None
+
+
+    _method = 'mixed'
+    @property
+    def method(self):
+        """
+
+            Method must be either 'mixed' or 'head'.
+
+            There are two different forms of Richards equation that differ
+            on how they deal with the non-linearity in the time-stepping term.
+
+            The most fundamental form, referred to as the
+            'mixed'-form of Richards Equation [Celia et al., 1990]
+
+            .. math::
+
+                \\frac{\partial \\theta(\psi)}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0
+                \quad \psi \in \Omega
+
+            where theta is water content, and psi is pressure head.
+            This formulation of Richards equation is called the
+            'mixed'-form because the equation is parameterized in psi
+            but the time-stepping is in terms of theta.
+
+            As noted in [Celia et al., 1990] the 'head'-based form of Richards
+            equation can be written in the continuous form as:
+
+            .. math::
+
+                \\frac{\partial \\theta}{\partial \psi}\\frac{\partial \psi}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0
+                \quad \psi \in \Omega
+
+            However, it can be shown that this does not conserve mass in the discrete formulation.
+
+
+        """
+        return self._method
+    @method.setter
+    def method(self, value):
+        assert value in ['mixed','head']
+        self._method = value
+
+    _doNewton = False
+    @property
+    def doNewton(self):
+        """Do a Newton iteration. If False, a Picard iteration will be completed."""
+        return self._doNewton
+    @doNewton.setter
+    def doNewton(self, value):
+        self._doNewton = value
+
+
+    def __init__(self, mesh, empirical):
+        Problem.__init__(self, mesh)
+        self.empirical = empirical
+
+
+    def getResidual(self, h, hn):
+        DIV  = self.mesh.faceDiv
+        GRAD = self.mesh.cellGrad
+        BC   = self.mesh.BC
+        AV   = self.mesh.AV
+        Dz   = self.mesh.Dz
+
+        T  = self.empirical.moistureContent(h)
+        dT = self.empirical.moistureContentDeriv(h)
+        Tn = self.empirical.moistureContent(hn)
+        K  = self.empirical.hydraulicConductivity(h)
+        dK = self.empirical.hydraulicConductivityDeriv(h)
+
+        aveK = 1./(AV*(1./K));
+
+        RHS = DIV*sdiag(aveK)*(GRAD*h+BC*P.bc) + Dz*aveK
+        if self.method is 'mixed':
+            r = (T-Tn)/dt - RHS
+        elif self.method is 'head':
+            r = dT*(h - hn)/dt - RHS
+
+        J = + dT/dt - DIV*sdiag(aveK)*GRAD
+        if self.doNewton:
+            DDharmAve = sdiag(aveK**2)*AV*sdiag(K**(-2)) * dK
+            J = J - DIV*sdiag(GRAD*h + BC*P.bc)*DDharmAve - Dz*DDharmAve
+
+        return r, J
+
+class Haverkamp(object):
+    """docstring for Haverkamp"""
+
+    empiricalModelName = "VanGenuchten"
+
+    theta_s = 0.430
+    theta_r = 0.078
+    alpha   = 0.036
+    beta    = 3.960
+    A       = 1.175e+06
+    gamma   = 4.74
+    Ks      = np.log(24.96)
+
+    def __init__(self, **kwargs):
+        setKwargs(self, **kwargs)
+
+    def moistureContent(self, h):
+        f = (self.alpha*(self.theta_s  -    self.theta_r  )/
+                        (self.alpha    + abs(h)**self.beta) + self.theta_r)
+        f[h > 0] = self.theta_s
+        return f
+
+    def moistureContentDeriv(self, h):
+        g = (self.alpha*((self.theta_s - self.theta_r)/
+             (self.alpha + abs(h)**self.beta)**2)
+             *(-self.beta*abs(h)**(self.beta-1)*np.sign(h)));
+        g[h >= 0] = 0
+        g = sdiag(g)
+        return g
+
+    def hydraulicConductivity(self, h):
+        f = np.exp(self.Ks)*self.A/(self.A+abs(h)**self.gamma)
+        if type(self.Ks) is np.ndarray and self.Ks.size > 1:
+            f[h >= 0] = np.exp(self.Ks[h >= 0])
+        else:
+            f[h >= 0] = np.exp(self.Ks)
+        return f
+
+    def hydraulicConductivityDeriv(self, h):
+        g = -(np.exp(self.Ks)*self.A*self.gamma*abs(h)**(self.gamma-1)*np.sign(h))/((self.A+abs(h)**self.gamma)**2)
+        g[h >= 0] = 0
+        g = sdiag(g)
+        #A
+        # dA = np.exp(self.Ks)/(self.A+abs(h)**self.gamma) - np.exp(self.Ks)*self.A/(self.A+abs(h)**self.gamma)**2;
+        #gamma
+        # dgamma = -(self.A*np.exp(self.Ks)*np.log(abs(h))*abs(h)**self.gamma)/(self.A + abs(h)**self.gamma)**2;
+        return g
+
+
+class VanGenuchten(object):
+    """
+
+    .. math::
+
+        \\theta(h) = \\frac{\\alpha (\\theta_s - \\theta_r)}{\\alpha + |h|^\\beta} + \\theta_r
+
+    Where parameters alpha, beta, gamma, A are constants in the media;
+    theta_r and theta_s are the residual and saturated moisture
+    contents; and K_s is the saturated hydraulic conductivity.
+
+    Celia1990
+
+    """
+
+    empiricalModelName = "VanGenuchten"
+
+    theta_s = 0.430
+    theta_r = 0.078
+    alpha   = 0.036
+    n       = 1.560
+    beta    = 3.960
+    I       = 0.500
+    Ks      = np.log(24.96)
+
+    def __init__(self, **kwargs):
+        setKwargs(self, **kwargs)
+
+    def moistureContent(self, h):
+        m = 1 - 1/self.n;
+        f = ((  self.theta_s  -  self.theta_r  )/
+             ((1+abs(self.alpha*h)**self.n)**m)   +  self.theta_r)
+        f[h > 0] = self.theta_s
+        return f
+
+    def moistureContentDeriv(self, h):
+        g = -self.alpha*self.n*abs(self.alpha*h)**(self.n - 1)*np.sign(self.alpha*h)*(1./self.n - 1)*(self.theta_r - self.theta_s)*(abs(self.alpha*h)**self.n + 1)**(1./self.n - 2)
+        g[h > 0] = 0
+        g = sdiag(g)
+        return g
+
+    def hydraulicConductivity(self, h):
+        alpha = self.alpha
+        I = self.I
+        n = self.n
+        Ks = self.Ks
+        m = 1 - 1/n
+
+        theta_e = 1/((1+abs(alpha*h)**n)**m)
+        f = np.exp(Ks)*theta_e**I* ( ( 1 - ( 1 - theta_e**(1/m) )**m )**2 )
+        if type(self.Ks) is np.ndarray and self.Ks.size > 1:
+            f[h >= 0] = np.exp(self.Ks[h >= 0])
+        else:
+            f[h >= 0] = np.exp(self.Ks)
+        return f
+
+    def hydraulicConductivityDeriv(self, h):
+        alpha = self.alpha
+        I = self.I
+        n = self.n
+        Ks = self.Ks
+        m = 1 - 1/n
+
+        g = I*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n))
+        g[h >= 0] = 0
+        g = sdiag(g)
+
+        #alpha
+        # dA = I*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n));
+        #n
+        # dn = 2*np.exp(Ks)*((np.log(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n))/n**2 + ((1/n - 1)*(((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/((1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1/n - 1))/(n**2*(1/n - 1)**2*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))))/(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2;
+        #I
+        # dI = np.exp(Ks)*np.log((abs(alpha*h)**n + 1)**(1/n - 1))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2;
+
+        return g
diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py
index 33c9a6b1..ec933f0a 100644
--- a/SimPEG/forward/__init__.py
+++ b/SimPEG/forward/__init__.py
@@ -2,3 +2,4 @@ from Problem import *
 import DCProblem
 from LinearProblem import LinearProblem
 import ModelTransforms
+import Richards
diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py
new file mode 100644
index 00000000..9c634e33
--- /dev/null
+++ b/SimPEG/tests/test_Richards.py
@@ -0,0 +1,66 @@
+import numpy as np
+import unittest
+from SimPEG.mesh import TensorMesh
+from TestUtils import OrderTest, checkDerivative
+from scipy.sparse.linalg import dsolve
+from SimPEG.forward import Richards
+
+
+class RichardsTests(unittest.TestCase):
+
+    def setUp(self):
+        pass
+        # a = np.array([1, 1, 1])
+        # b = np.array([1, 2])
+        # c = np.array([1, 4])
+        # self.mesh2 = TensorMesh([a, b], np.array([3, 5]))
+        # self.mesh3 = TensorMesh([a, b, c])
+
+    def test_VanGenuchten_moistureContent(self):
+        vanG = Richards.VanGenuchten()
+        def wrapper(x):
+            return vanG.moistureContent(x), vanG.moistureContentDeriv(x)
+        passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
+        self.assertTrue(passed,True)
+
+    def test_VanGenuchten_hydraulicConductivity(self):
+        hav = Richards.VanGenuchten()
+        def wrapper(x):
+            return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
+        passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
+        self.assertTrue(passed,True)
+
+    def test_VanGenuchten_hydraulicConductivity_FullKs(self):
+        n = 50
+        hav = Richards.VanGenuchten(Ks=np.random.rand(n))
+        def wrapper(x):
+            return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
+        passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False)
+        self.assertTrue(passed,True)
+
+    def test_Haverkamp_moistureContent(self):
+        hav = Richards.Haverkamp()
+        def wrapper(x):
+            return hav.moistureContent(x), hav.moistureContentDeriv(x)
+        passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
+        self.assertTrue(passed,True)
+
+    def test_Haverkamp_hydraulicConductivity(self):
+        hav = Richards.Haverkamp()
+        def wrapper(x):
+            return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
+        passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
+        self.assertTrue(passed,True)
+
+    def test_Haverkamp_hydraulicConductivity_FullKs(self):
+        n = 50
+        hav = Richards.Haverkamp(Ks=np.random.rand(n))
+        def wrapper(x):
+            return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
+        passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False)
+        self.assertTrue(passed,True)
+
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst
index 068e30da..7fea7b5c 100644
--- a/docs/api_Problem.rst
+++ b/docs/api_Problem.rst
@@ -26,3 +26,10 @@ Linear Problem
     :members:
     :undoc-members:
 
+Richards Problem
+****************
+
+.. automodule:: SimPEG.forward.Richards
+    :members:
+    :undoc-members:
+