  Simulation of geothermal systems is challenging due to coupled physical
processes in highly heterogeneous media. Combining the exponential
Rosenbrock--Euler and Rosenbrock-type methods with control-volume (two-point
flux approximation) space discretizations leads to efficient numerical
techniques for simulating geothermal systems. In terms of efficiency and
accuracy, the exponential Rosenbrock--Euler time integrator has advantages over
standard time-dicretization schemes, which suffer from time-step restrictions
or excessive numerical diffusion when advection processes are dominating. Based
on linearization of the equation at each time step, we make use of matrix
exponentials of the Jacobian from the spatial discretization, which provide the
exact solution in time for the linearized equations. This is at the expense of
computing the matrix exponentials of the stiff Jacobian matrix, together with
propagating a linearized system. However, using a Krylov subspace or Leja
points techniques make these computations efficient. The Rosenbrock-type
methods use the appropriate rational functions of the Jacobian from the spatial
discretization. The parameters in these schemes are found in consistency with
the required order of convergence in time. As a result, these schemes are
A-stable and only a few linear systems are solved at each time step.
  The efficiency of the methods compared to standard time-discretization
techniques are demonstrated in numerical examples.
