  The bottomonium spectrum is computed in dynamical 2+1 flavor lattice QCD,
using NRQCD for the b quarks. The main calculations in this work are based on
gauge field ensembles generated by the RBC and UKQCD collaborations with the
Iwasaki action for the gluons and a domain-wall action for the sea quarks.
Lattice spacing values of approximately 0.08 fm and 0.11 fm are used, and
simultaneous chiral extrapolations to the physical pion mass are performed. As
a test for gluon discretization errors, the calculations are repeated on two
ensembles generated by the MILC collaboration with the Luscher-Weisz gauge
action. Gluon discretization errors are also studied in a lattice potential
model using perturbation theory for four different gauge actions. The
nonperturbative lattice QCD results for the radial and orbital bottomonium
energy splittings obtained from the RBC/UKQCD ensembles are found to be in
excellent agreement with experiment. To get accurate results for spin
splittings, the spin-dependent order-v^6 terms are included in the NRQCD
action, and suitable ratios are calculated such that most of the unknown
radiative corrections cancel. The cancellation of radiative corrections is
verified explicitly by repeating the calculations with different values of the
couplings in the NRQCD action. Using the lattice ratios of the S-wave hyperfine
and the 1P tensor splitting, and the experimental result for the 1P tensor
splitting, the 1S hyperfine splitting is found to be
60.3+-5.5(stat)+-5.0(syst)+-2.1(exp) MeV, and the 2S hyperfine splitting is
predicted to be 23.5+-4.1(stat)+-2.1(syst)+-0.8(exp) MeV.
