  We use $n_f=2+1$ Wilson-clover gauge-field ensembles from the CLS consortium
in a Lattice NRQCD setup to predict the binding energy of a $I(J^P)=0(1^+)$
$ud\bar{b}\bar{b}$ tetraquark and a $\frac{1}{2}(1^+)$ $\ell s\bar{b}\bar{b}$
tetraquark. We determine the binding energies with respect to the relevant
$BB^*$ and $B_sB^*$ thresholds respectively to be $112.0(13.2)$ MeV for the
$ud\bar{b}\bar{b}$, and $46.4(12.3)$ MeV for the $\ell s\bar{b}\bar{b}$. We
also determine the ground-state $J^P=0^+$ $B_{s0}^*$ and $1^+$ $B_{s1}$ mesons
to lie $75.4(14.0)$ and $78.7(13.9)$ MeV below the $BK$ and $B^*K$ thresholds
respectively. Our errors are entirely dominated by systematics due to
discretisation effects. To achieve these measurements, we performed a neural
network based nonperturbative tuning of the Lattice NRQCD Hamiltonian's
parameters against the basic bottomonium spectrum. For all lattice spacings
considered we can reproduce the continuum splittings of low-lying bottomonia.
It is worth remarking that our nonperturbative tuning parameters deviate from 1
by significant amounts, particularly the term $c_2$.
