def conduction(z, T1, qb, lam, z0):
    """
    Compute conduction temperature field from Domenico (1978).

    Formula (11).

    Domenico, P. A., & Palciauskas, V. V., Theoretical analysis of forced
    convective heat transfer in regional ground-water flow,
    Geological Society of America Bulletin, 84(12), 3803 (1973).
    http://dx.doi.org/10.1130/0016-7606(1973)84<3803:taofch>2.0.co;2

    Parameters
    ----------
    z : Numpy array
        Three dimensional array of z-coordinates.

    T1 : float
        background temperature

    qb : float
        base heat flow

    k : float
        permeability

    lam : float
        thermal conductivity of the mixture
        of matrix and fluid

    z0 : float
        Parameter in formula. Characteristic length
        in z-direction. Upper surface of flow.

    Returns
    ----------
    T: Numpy array
        3D array, temperature field.
    """

    # Temperature gradient
    dTdz = qb / lam

    # Temperature
    T = (
        T1 +  # Constant background
        dTdz * (z0 - z)  # Simple gradient
    )

    return T