  The clustered nature of star formation should produce a high degree of
structure in the combined phase and chemical space in the Galactic disk. To
date, observed structure of this kind has been mostly limited to bound clusters
and moving groups. In this paper we present a new dynamical model of the
Galactic disk that takes into account the clustered nature of star formation.
This model predicts that the combined phase and chemical space is rich in
substructure, and that this structure is sensitive to both the precise nature
of clustered star formation and the large-scale properties of the Galaxy. The
model self-consistently evolves 4 billion stars over the last 5 Gyr in a
realistic potential that includes an axisymmetric component, a bar, spiral
arms, and giant molecular clouds (GMCs). All stars are born in clusters with an
observationally-motivated range of initial conditions. As direct
\textit{N}-body calculations for billions of stars is computationally
infeasible, we have developed a method of initializing star cluster particles
to mimic the effects of direct \textit{N}-body effects, while the actual orbit
integrations are treated as test particles within the analytic potential. We
demonstrate that the combination of chemical and phase space information is
much more effective at identifying truly co-natal populations than either
chemical or phase space alone. Furthermore, we show that co-moving pairs of
stars are very likely to be co-natal if their velocity separation is $< 2$ km
s$^{-1}$ and their metallicity separation is $< 0.05$ dex. The results
presented here bode well for harnessing the synergies between \textit{Gaia} and
spectroscopic surveys to reveal the assembly history of the Galactic disk.
