# prefflow_snowpack.py

This is a code imitating the preferential flow scheme developed in Wever et al. (2016) Water flow in firn with dual domain approach: - Richards Equation in Matrix Flow domain and Preferential Flow domain

Few differences with Wever 2016 as: -Use of a constant bigF value (part of the pore space allocated to each domain) -Runoff function from Zuo and Oerlemans 1996 -Possible to apply PFfreezing if cold wave penetrates from the surface -Possible to build up aquifer at end of domain -Use of upstream weighted mean to determine hydraulic conductivity at interfaces avoids oscillations for large mesh size -We use a changing bottom boundary if saturated layers accumulate: don’t solve RE for the saturated layers at end of the firn column - don’t solve RE in dry part of the domain

Density of last layer of the domain should always be >= 830