  We cast the classical Yang-Baxter equation (CYBE) in a variational context
for the first time, by relating it to the theory of Lagrangian multiforms, a
framework designed to capture integrability in a variational fashion. This
provides a significant connection between Lagrangian multiforms and the CYBE,
one of the most fundamental concepts of integrable systems. This is achieved by
introducing a generating Lagrangian multiform which depends on a skew-symmetric
classical $r$-matrix with spectral parameters. The multiform Euler-Lagrange
equations produce a generating Lax equation which yields a generating zero
curvature equation. The CYBE plays a role at three levels: 1) It ensures the
commutativity of the flows of the generating Lax equation; 2) It ensures that
the generating zero curvature equation holds; 3) It implies the closure
relation for the generating Lagrangian multiform. The specification of an
integrable hierarchy is achieved by fixing certain data: a finite set $S\in
CP^1$, a Lie algebra $\mathfrak{g}$, a $\mathfrak{g}$-valued rational function
with poles in $S$ and an $r$-matrix. We show how our framework is able to
generate a large class of ultralocal integrable hierarchies by providing
several known and new examples pertaining to the rational or trigonometric
class. These include the Ablowitz-Kaup-Newell-Segur hierarchy, the sine-Gordon
(sG) hierarchy and various hierachies related to Zakharov-Mikhailov type models
which contain the Faddeev-Reshetikhin (FR) model and recently introduced
deformed sigma/Gross-Neveu models as particular cases. The versatility of our
method is illustrated by showing how to couple integrable hierarchies together
to create new examples of integrable field theories and their hierarchies. We
provide two examples: the coupling of the nonlinear Schr\"odinger system to the
FR model and the coupling of sG with the anisotropic FR model.
