We formally derive interface conditions for modeling thin inclusions in heterogeneous diffusion problems expressed in the form of the divergence of a flux. By integrating the governing equations within the inclusions, we establish that the resulting interface conditions are of Wentzell-type for the flux jump and of Robin-type for the flux average. The condition on the flux jump condition involves an unconventional tangential diffusion operator applied to the average solution across the interface. The corresponding weak formulation is introduced, providing a framework that is readily applicable to finite element discretizations. Extensive numerical validation does not only demonstrate the robustness of the proposed modeling technique but also, and in particular, its ability to accomodate complex networks of thin inclusions without any need for an explicit geometric representation of the latter. Our tests also demonstrate that this novel approach is able to handle extremely strong differences in the material properties between the inclusions and their embedding background as they are indeed typical for flow problems involving fractured media.