  We classify and construct models for two-dimensional (2D) interacting
fermionic symmetry-protected topological (FSPT) phases with general finite
Abelian unitary symmetry $G_f$. To obtain the classification, we couple the
FSPT system to a dynamical discrete gauge field with gauge group $G_f$ and
study braiding statistics in the resulting gauge theory. Under reasonable
assumptions, the braiding statistics data allows us to infer a potentially
complete classification of 2D FSPT phases with Abelian symmetry. The FSPT
models that we construct are simple stacks of the following two kinds of
existing models: (i) free-fermion models and (ii) models obtained through
embedding of bosonic symmetry-protected topological (BSPT) phases.
Interestingly, using these two kinds of models, we are able to realize almost
all FSPT phases in our classification, except for one class. We argue that this
exceptional class of FSPT phases can never be realized through models (i) and
(ii), and therefore can be thought of as intrinsically interacting and
intrinsically fermionic. The simplest example of this class is associated with
$\mathbb Z_4^f\times\mathbb Z_4\times\mathbb Z_4$ symmetry. We show that all 2D
FSPT phases with a finite Abelian symmetry of the form $ \mathbb Z_2^f\times G$
can be realized through the above models (i), or (ii), or a simple stack of
them. Finally, we study the stability of BSPT phases when they are embedded
into fermionic systems.
