  The valleys in hexagonal two-dimensional systems with broken inversion
symmetry carry an intrinsic orbital magnetic moment. Despite this, such systems
possess zero net magnetization unless additional symmetries are broken, since
the contributions from both valleys cancel. A nonzero net magnetization can be
induced through applying both uniaxial strain to break the rotational symmetry
of the lattice and an in-plane electric field to break time-reversal symmetry
owing to the resulting current. This creates a magnetoelectric effect whose
strength is characterized by a magnetoelectric susceptibility, which describes
the induced magnetization per unit applied in-plane electric field. Here, we
predict the strength of this magnetoelectric susceptibility for Bernal-stacked
bilayer graphene as a function of the magnitude and direction of strain, the
chemical potential, and the interlayer electric field. We estimate that an
orbital magnetization of ~5400 $\mu_{\text{B}}/\mu\text{m}^2$ can be achieved
for 1% uniaxial strain and a 10 $\mu\text{A}$ bias current, which is almost
three orders of magnitude larger than previously probed experimentally in
strained monolayer MoS$_2$. We also identify regimes in which the
magnetoelectric susceptibility not only switches sign upon reversal of the
interlayer electric field but also in response to small changes in the carrier
density. Taking advantage of this reversibility, we further show that it is
experimentally feasible to probe the effect using scanning magnetometry.
