Constitutive laws relating fluid potentials and fluxes in a nonlinear manner are common in several porous media applications, including biological and reactive flows, poromechanics, and fracture deformation. Compared to the standard, linear Darcy’s law, such enhanced flux relations increase both the degree of nonlinearity, and, in the case of multiphysics simulations, coupling strength between processes. While incorporating the nonlinearities into simulation models is thus paramount for computational efficiency, correct linearization, as is needed for incorporation in Newton’s method, is challenging from a practical perspective. The standard approach is therefore to ignore nonlinearities in the permeability during linearization. For finite volume methods, which are popular in porous media applications, complete linearization is feasible only for the simplest flux discretization, namely the two-point flux approximation. We introduce an approximated linearization scheme for finite volume methods that is exact for the two-point scheme and can be applied to more advanced and accurate discretizations, exemplified herein by a multi-point flux stencil. We test the new method for both nonlinear porous media flow and several multiphysics simulations. Our results show that the new linearization consistently outperforms the standard approach. Moreover our scheme achieves asymptotic second order convergence of the Newton iterations, in contrast to the linear convergence obtained with the standard approach.