Accurate representation and estimation of seismic sources are essential to the seismic inversion problem. General sources can be approximated by a truncated series of multipoles depending on the source anisotropy. Most research in the joint inversion of source and medium parameters assumes seismic sources can be modeled as isotropic point-sources resulting in an inability to fit the anisotropy observed in data, ultimately impacting the recovery of medium parameters. In this thesis I lay the groundwork for joint source-medium parameter inversion with potentially anisotropic seismic sources via full waveform inversion through three key contributions: a mathematical and computational framework for the modeling and inversion of sources via multipoles, construction and analysis of discretizations of multipole sources on regular grids, and preconditioners based on fractional time derivative/integral operators for the ill-conditioned source estimation subproblem. As an application of my multipole framework, I also study the efficacy of multipoles in modeling the airgun array source, the most common type of active source in marine seismic surveying. Inversion results recovered a dominating isotropic component of the multipole source model that accounted for 84% of the observed radiation pattern. An extra 10% of the observed output pressure field can be explained when incorporating dipole terms in the source representation, thus motivating the use of multipoles to capture source anisotropy.