High-fidelity physics-based numerical models have often been utilized for subsurface modeling in the oil and gas industry. The thermo-hydro-chemo-mechanical interactions in the subsurface processes are incorporated by means of coupled-field balance laws in these models. The resulting equations are discretized and then solved on a computer using a fit-for-purpose numerical method. However, all numerical methods require that the subsurface system be appropriately meshed into non-overlapping elements. Unfortunately, generating such meshes for subsurface systems is a non-trivial task. This is because most subsurface systems are geometrically complex and exhibit a high degree of heterogeneity. In recent years, a novel class of approaches known as Physics Informed Neural Networks (PINNs) is getting increasingly popular. In PINNs, the governing equations are solved by casting the problem into a PDE-constrained optimization problem. The loss functionals are formed from the residuals of physics equations and boundary and/or initial conditions. In contrast to purely data-driven modelling approaches, these methods are informed by the underlying physics and, therefore, require lesser training data. Barring a few exceptions, these methods have yet to be thoroughly studied in the context of heterogeneous systems. This study investigates the potential of Physics Informed Neural Networks (PINNs) in modeling heterogeneous systems and proposes a novel architecture of PINNs, referred to as "Interface-PINNs," that is specifically designed to address the interfacial discontinuities within the heterogeneous domain. Additionally, the study explores another framework of PINNs, known as Variational PINNs (V-PINNs), that minimizes the potential energy of a system instead of the governing differential equations in the strong form. A comprehensive comparative analysis between the two frameworks (PINNs and V-PINNs) is conducted. Several benchmark numerical examples in 1-D and 2-D have been studied to investigate the efficacy of both PINNs and V-PINNs with the proposed "Interface- PINNs" architecture in capturing strong and weak continuities at interfaces. The results of the study demonstrate that both PINNs and V-PINNs, coupled with the proposed "Interface-PINNs" architecture, are capable of accurately capturing discontinuities at the interfaces, thereby producing precise approximations of the closed-form solutions.