The first releases of porousMultiphaseFoam proposed an open-source software suite to solve the equations for multiphase flow (generalized Darcy’s law) in porous media or groundwater flows (Richards' equation) by taking advantage of OpenFOAM, a finite volume platform with automatic discretization on three-dimensional unstructured grids and good parallel efficiency. Recently, the porousMultiphaseFoam toolbox has been confronted with complex cases of fast water flows and solute transfers in realistic hydrological configurations with variable forcing conditions (heterogeneous infiltration and local tracer injection). Several developments have been carried out to make it possible to simulate those cases, which extend the toolbox with: (i) a set of solvers dedicated to groundwater flows, including coupled water flow and solute transport and simplified 2D approaches, (ii) improved numerical techniques for problems with strong non-linearities, (iii) libraries/executables for pre-processing of input data (geographical information and time-variable forcing terms) and (iv) passive or coupled scalar transport (tracer) with groundwater solvers that support any number of species. New solvers are validated on several (un-)saturated configurations by a direct comparison with a well validated finite element code. Program summaryProgram title: porousMultiphaseFoam v2107CPC Library link to program files:https://doi.org/10.17632/hphn58ksfy.1Developer's repository link:https://github.com/phorgue/porousMultiphaseFoam.gitLicensing provisions: GNU GPLv3Programming language: C++ / OpenFOAMJournal reference of previous version: Computer Physics Communications 187 (2015) 217-226Does the new version supersede the previous version?: YesReasons for the new version: Multiple new features of the toolbox dedicated to groundwater flows (with pre-processing tools)Summary of revisions:–Solvers:–groundwaterFoam/stationaryGroundwaterFoam: new unsaturated flow (transient/stationary) solvers in porous media (Richards' equation) using Picard or Newton's method for linearization.–groundwater2DFoam/stationaryGroundwater2DFoam: new 2D (transient/stationary) solver for saturated flow in porous media.–porousScalarTransportFoam/porousScalarTransport2DFoam: new passive multi-scalar transport solvers on pre-computed flow fields (3D unsaturated or 2D saturated fields).–groundwaterTransportFoam/groundwaterTransport2DFoam: new solvers coupling flow modeling (2D saturated or 3D unsaturated) with multi-scalar tracers transports.–Improved numerical techniques–Timestep management: time-stepping is now based on the temporal truncation of the time scheme used (backward Euler 1st and 2nd order, Crank-Nicolson).–Linearization techniques: problems of unsaturated flow with strong non-linearities can now be solved using an algorithm based on Newton's method.–Libraries:–numericalMethods: the EulerD3dt3Scheme class allows computation of the 3rd-order time derivative to estimate the truncation error and the JacobianMatrix class computes the Jacobian using finite differencing (expensive method used only for debugging).–toolsGIS: new dynamic library for managing GIS (Geographical Information System) files and event files. Events are variable and localized source terms for fluid or scalar transport.–Ippisch: an improved two-phase model for relative permeability and capillary pressure.–porousBoundaryConditions: several new boundary conditions have been added to handle heterogeneous pressure head boundary conditions and variable fluid or tracer injections (using event files).–Utilities:–setBoundaryHeadPressure: pressure head initializer for boundary conditions using uniform value of GIS files.–setFieldsFromDEM/setFieldsFromXY: heterogeneous initializer for internal fields, similar to the classical OpenFOAM setFields, but using ordered (DEM) or unordered (XY) grid files.–darcyFoam: simple utility to calculate a single-phase flow in porous media.Nature of problem: The porousMultiphaseFoam toolbox was initially released with generic two-phase solvers for flow in porous media. Developments were carried out to model more specifically hydrogeological flows by including a groundwaterFoam solver for Richards' equation. However, many features required to simulate realistic configurations, such as event management for variable and heterogeneous forcing terms (infiltration / scalar source term) or suitable numerical techniques for complex configurations (high rate injection for unsaturated flow for example), were still missing. Moreover, the 2D approach with integration of vertical properties (for flow and tracer) was not present even though this is an indispensable tool for groundwater flow modeling due to its ease of implementation and its low computation cost.Solution method: Improved numerical techniques have been developed to solve realistic flows by handling strong non-linearities with Newton's method, and a wide range of time scales with time-stepping based on time truncation errors. Various solvers for groundwater flows have been developed to handle cases under different assumptions, including the aforementioned 2D approaches. Variants of these solvers were also developed to simultaneously solve the passive transport of any number of scalars (tracers) by the flow. A dynamic library has been developed to pre-process input data: event information for variable conditions, GIS file for spatial information.Additional comments including restrictions and unusual features: porousMultiphaseFoam v2107 requires an installation of OpenFOAM-v9. Dedicated branches are available on the developer's repository for compatibility with other versions of OpenFOAM. The following versions are compatible with the current version of porousMultiphaseFoam: v7, v8, v1906, v2006 and v2106. Older branches may be found in the repository link but do not contain all features.