...the quantity
$$\vec{v}(\vec{k})=\frac{1}{\hbar}\nabla_{\vec{k}}E_n(\vec{k}),$$ ...
...
$$= \frac{1}{m}\langle\psi_{n,\vec{k}}|-i\hbar\nabla|\psi_{n,\vec{k}}\rangle.$$
Yes, this is derived in Ashcroft and Mermin's Appendix E.
So far there are no wavepackets and electrons are described by delocalized Bloch waves.
Yes, and remember that you can roughly think of Bloch waves in a periodic potential as analogous to plane waves for free particles (no potential).
However, when one talks about the semiclassical dynamics of electrons, one describes the electrons not in terms of delocalized Bloch waves but by a somewhat localized wavepacket. Why is this so?
Both Bloch waves and plane waves are completely delocalized; they have non-zero amplitudes everywhere in space, so they can not describe a localized particle.
In order to describe a localized particle in a periodic potential one can use a wave packet of Bloch waves, just like one can use a wave packet of plane waves to describe a localized free particle. Both Bloch waves and plane waves are complete sets of functions, so the choice is a matter of convenience. Bloch waves can be more convenient to use when the problem statement includes a periodic (non-constant) potential.
In short, the question is: Why wavepackets?
Because we have to use wave packets to describe localized particles. Plane waves and Bloch waves are delocalized and unphysical, in the sense that they are unnormalizable (unless you put them in a box). For example, with plane waves, the probability amplitude to find the particle anywhere in space is:
$$
P_{pl}(\vec r) = |e^{i\vec k \cdot \vec r}|^2 = 1\;.
$$
Or, if we use "box normalization" where the box has a size $V$, the probability amplitude is still constant and the same everywhere in space:
$$
P_{pl,V}(\vec r) = \frac{1}{V}\;.
$$
Similarly, with Bloch waves, the probability amplitude is:
$$
P_{bl}(\vec r) = |e^{i\vec k \cdot \vec r}u(\vec r)|^2 = |u(\vec r)|^2\;,
$$
which has the same non-zero values in any unit cell in space, since $u(\vec r +\vec R)=u(\vec r)$.
We have grown used to thinking of plane waves as describing free particles, but really they do not; they happen to solve the free-particle time-independent Schrödinger equation, but they are not normalizable (unless forced into a box). But wave packets made up of plane waves can be normalized. Similarly, wave packets made up of Bloch waves can be normalized.
Then, why $\vec{v}=(1/\hbar)\nabla_{\vec{k}}E_n(\vec{k})$ is its group velocity of that packet?
Suppose you have some physical/normalizable wave packet $\Psi(\vec r, 0)$ where
$$
\Psi(\vec r, 0 ) = \sum_{n\vec k}c_{n\vec k}e^{i\vec k \cdot \vec r}u_{n\vec k }(\vec r)\;.
$$
The time evolution is straight forward to figure out since we expanded in terms of Bloch waves:
$$
\Psi(\vec r, t) = \sum_{n\vec k}c_{n\vec k}e^{-iE_{n}(\vec k) t}e^{i\vec k \cdot \vec r}u_{n,\vec k }(\vec r)\;.
$$
Now suppose that the packet is peaked about some specific $\vec k = \vec k_0$ and $n=n_0$ and write:
$$
\Psi(\vec r, t) \approx e^{-i(E_{n_0}(\vec k_0)+\left.\frac{\partial E_{n_0}}{\partial \vec k}\right|_{\vec k_0}\cdot (\vec k - \vec k_0)+\ldots) t}e^{i\vec k_0 \cdot \vec r}\sum_{\vec k} c_{n_0,\vec k}e^{i(\vec k - \vec k_0)\cdot r}u_{n_0,\vec k }(\vec r)\;,
$$
where we can think of $c_{n_0\vec k}$ as something like $e^{-\alpha(\vec k - \vec k_0)^2}$ so that only $\vec k$ values near $\vec k_0$ are important for large $\alpha$. Continuing, we can thus write
$$
\Psi(\vec r, t)\approx e^{-i\vec k_0\cdot (\vec r + \vec v_0 t)}e^{-i E_{n_0}(\vec k_0)t}A(\vec r)\;,
$$
where
$$
\vec v_0 \equiv \left.\frac{\partial E_{n_0}(\vec k)}{\partial \vec k}\right|_{\vec k_0}
$$
and where $A(\vec r)$ is some function that is localized and we can think of it as looking something like
$$
A(\vec r) \approx e^{-\frac{r^2}{2\alpha}}\;.
$$
So, our wave packet "looks like" a plane wave, but with $\vec r$ shifted to $\vec r + \vec v_0 t$ and some overall envelope function $A(\vec r)$ keeping it localized.
