$
\newcommand\q\mathbf
\newcommand\grade[1]{\langle#1\rangle}
\newcommand\PD[2]{\frac{\partial#1}{\partial#2}}
\newcommand\re{\mathrm{re}}
$
Edit:
I realized that my original post likely didn't address your actual question: how to compute $\PD{}{q_a}\q Q_1^2$, where I assume you mean $a=1,2,3$. With partial derivatives, we always have to keep in mind what is varying and what is being held constant. In this case, since $(x, y, z)$ is one set of variables and $(q_1, q_2, q_3)$ is another, I would assume in taking a derivative like $\PD{}{q_1}\q Q_1^2$ that $q_2, q_3$ are held constant and $x, y, z$ are allowed to vary as function of $q_1, q_2, q_3$. So then, since $\q Q_1 = q_1\q i + q_2\q j + q_3\q k$, we see
$$
\PD{\q Q_1}{q_1} = \PD{}{q_1}(q_1\q i + q_2\q j + q_3\q k) = \q i
$$
and so we get
$$
\PD{}{q_1}\q Q_1^2 = \PD{\q Q_1}{q_1}\q Q_1 + \q Q_1\PD{\q Q_1}{q_1} = \q i\q Q_1 + \q Q_1\q i = \frac12\re[\q i\q Q_1],
$$
the last equality following because both $\q i$ and $\q Q_1$ are pure imaginary. In the case that you want to compute something like $\PD{}x\q Q_1^2$, first we note that
$$
q_1 = \frac x{2\sqrt z},\quad q_2 = \frac y{2\sqrt z},\quad q_3 = \frac12\sqrt z,
$$
and then by the chain rule we get
$$
\PD{}x\q Q_1^2
= \sum_{a=1}^3\PD{q_a}x\PD{}{q_a}\q Q_1^2
= \frac1{2\sqrt z}\PD{}{q_1}\q Q_1^2
= \frac1{4\sqrt z}\re[i\q Q_1].
$$
Addendum to my original post:
Having spent some time looking at your second link, I think I understand what is going on now. Firstly, the operator $\nabla_{\q Q_1}$ from your first link is not the gradient $\nabla$ from your second link. It can't be, since in your second link the gradient is a 4D vector of quaternions (i.e. an element of $\mathbb H^4$) and the identity $\nabla_{\q Q_1}\q Q_1^2 = 2\q Q_1$ wouldn't make any sense. What the operator $\nabla_{\q Q_1}$ appears to correspond to is the $\PD{}q$ from your second link:
$$
\nabla_{\q Q_1} = \PD{}q = \frac14\left[q_0 - \q i\PD{}{q_1} - \q j\PD{}{q_2} - \q k\PD{}{q_3}\right],
$$
where I've used $0,1,2,3$ instead of their $a,b,c,d$. This is exactly $\frac14$ times the even multivector derivative of $\mathrm{Cl}_3(\mathbb R)$ I refer to below. However, as I also pointed out below, in this case the identity $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = A$ is incorrect, and should instead be $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = A^*$. In fact, the only consistent possibilities are (for imaginary $\q Q_1$)
$$
\nabla_{\q Q_1}\q Q_1^2 = 2\q Q_1,\quad \nabla_{\q Q_1}\re[\q Q_1^*\q A] = A^*
$$
or
$$
\nabla_{\q Q_1}\q Q_1^2 = 2\q Q_1^*,\quad \nabla_{\q Q_1}\re[\q Q_1^*\q A] = A
$$
in light of the fact that $\q Q_1^2 = \re[\q Q_1\q Q_1]$. (And what these correspond to is taking $\re[\q A\q B]$ or $\re[\q A^*\q B]$ as our desired inner product.)
On a tangent, I think care needs to be taken with the $\PD{}q, \PD{}{q^i},\PD{}{q^j},\PD{}{q^k}$ from your second link. While it is indeed true that, for example, $\PD q{q^i} = \PD{q^i}q = 0$, they are not independent, for consider
$$
\PD{}qiq^i = -\PD{}qi^2qi = \PD{}qqi = i.
$$
As far as I can tell, the authors do not comment on this.
Original post:
I'm having a hard time figuring out how $\nabla_{\q Q_1}$ is being defined, but if it's reasonable enough then this should be on the right track.
We need only consider arbitrary pure imaginary $\q Q_1$. The trick is that $\q Q_1^2 = \re[\q Q_1^2]$ and $\re[\q A\q B] = \re[\q B\q A]$. Hence, using the product rule and dots to indicate what variables are being differentiated,
$$
\nabla_{\q Q_1}\q Q_1^2
= \nabla_{\q Q_1}\re[\q Q_1^2]
= \dot\nabla_{\q Q_1}\re[\dot{\q Q_1}\q Q_1] + \dot\nabla_{\q Q_1}\re[\q Q_1\dot{\q Q_1}]
= 2\dot\nabla_{\q Q_1}\re[\dot{\q Q_1}\q Q_1].
$$
However, the second identity given in your first link, $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = \q A$, would then suggest that
$$
\nabla_{\q Q_1}\q Q_1^2 = 2\dot\nabla_{\q Q_1}\re[\dot{\q Q_1}\q Q_1] = 2\q Q_1^* = -2\q Q_1
$$
This appears to be an inconsistency. Again, on one hand I might not be understanding how $\nabla_{\q Q_1}$ is defined, but on the other hand I didn't assume very much. If we define $\nabla_{\q Q_1}$ for $\q Q_1 = q_0 + q_1\q i + q_2\q j + q_3\q k$ as
$$
\nabla_{\q Q_1} = \PD{}{q_0} - \q i\PD{}{q_1} - \q j\PD{}{q_2} - \q k\PD{}{q_3},
$$
then I am certain that my derivation is correct up until application of the $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = \q A$ identity. In fact, in this case I would say that $\nabla_{\q Q_1}\q Q_1^2 = 2\q Q_1$ is correct and $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = \q A$ is wrong.
I should point out that $\re[\q Q_1^*\q A]$ is exactly the Euclidean inner product on $\q Q_1$ and $\q A$ considered as 4D vectors, and the two identities $\nabla_{\q Q_1}\q Q_1^2 = 2\q Q_1$ and $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = \q A$ are identities of the vector derivative in geometric calculus, calculus done over real Clifford algebras (but the $\q Q_1^2$ identity would require using the Clifford product instead of the quaternion product). See Chapter 2 of Clifford Algebra to Geometric Calculus (1984) by David Hestenes and Garret Sobczyk, or Chapter 6 of Geometric Algebra for Physicists (2003) by Chris Doran and Anthony Lasenby; I recommend the latter as an introduction.
What does work with quaternion multiplication in this framework, though, is the case $\nabla_{\q Q_1}\re[\q Q_1^*\q A] = \q A^*$. This corresponds to treating quaternions as elements of the even subalgebra of $\mathrm{Cl}_3(\mathbb R)$ and using the even multivector derivative; this is arguably the more natural approach. With this approach, for pure imaginary $\q Q_1$ we do get the requisite $\nabla_{\q Q_1}\q Q_1^2 = 2\q Q_1$. In this case we can express the derivative of an arbitrary $\q Q_1$ as
$$
\nabla_{\q Q_1} = \PD{}{q_0} - \q i\PD{}{q_1} - \q j\PD{}{q_2} - \q k\PD{}{q_3}.
$$
