Electromagnetic metasurfaces are 2D arrays of sub-wavelength resonating particles whose microscopic properties can be tailored to achieve a desired macroscopic field scattering response. Determining the scattered fields from large metasurfaces consisting of small resonators constitutes a multi-scale problem. To remedy this, metasurfaces are modeled as zero thickness sheets described using surface susceptibilities in conjunction with the GSTCs. Here, the GSTCs are implemented in an FMM-IE solver to simulate metasurfaces embedded in complex environments. While the FMM-IE-GSTC solver works well, the standard dipolar surface susceptibility model of metasurfaces is inadequate for modeling structures exhibiting spatial dispersion. These spatially dispersive structures can be modeled using surface susceptibilities that are rational polynomial functions of the transverse wave-vector leading to the extended GSTCs. In this thesis, the extended GSTCs are further developed to model non-uniform metasurfaces providing a computationally efficient zero thickness model for practical metasurfaces which is then integrated into the IE-GSTC solver.