write test for global mass conservation
We have to test if global water mass is conserved. To do so, we need a new local operator which computes V_w = \sum_{T} \int_T \theta_w
. We can then compare the temporal change of this quantity to the fluxes into and out of the domain, d V_w / dt = \sum_{F} \int_F j_N
.