This work presents a number of methodologies for determining the scattered Electromagnetic (EM) fields from metasurfaces at various simulation scales. The primary motivation for this work is the need to create a suite of simulation techniques that enable the multi-scale simulation of large scale EM propagation problems involving metasurfaces. For all of the methodologies, the scattering behavior of metasurfaces are characterized by the surface susceptibility tensors that appear in the generalized sheet transition conditions (GSTCs), which gives a physically motivated description of the surface. Using the GSTCs, these shared constitutive parameters are integrated into several different modelling techniques, each with their own advantages, allowing for the simulation of metasurfaces at several simulation scales. This work will present several simulation frameworks that have been developed, outlining their various advantages and limitations. This includes: a modification of Finite Difference Time Domain (FDTD) method to integrate non-paraxial scattering behavior of metasurfaces to a standard Yee cell grid; the integration of metasurfaces as surface boundary conditions in the Boundary Element Method (BEM); the development of ray diffraction methods to calculate the scattered fields from metasurfaces characterized by local periodic parameters, and; the expansion of the ray diffraction methods to account for surface diffraction off of curved topographies. For each method, a few representative examples are included to verify their accuracy. The methods outlined above are demonstrated to be self-consistent with each other and with other techniques available in the present literature and offer a wide range of coverage to simulate metasurfaces accurately for any given simulation scale.