My answer to (1) is quick but some brief bits of background first.
The orbit-stabilizer theorem says that when a nice topological group has a nice continuous action on a topological space, any choice of basepoint yields a bundle $\mathrm{Stab}(x)\to G\to\mathrm{Orb}(x)$. The fibration is $\pi(g):=gx$ with fibers $\pi^{-1}(gx)=g\mathrm{Stab}(x)$; that is, the fibers are cosets of stabilizers. Also helpful to keep in mind $\mathrm{Stab}(gx)=g\mathrm{Stab}(x)g^{-1}$, so points in the same orbit have conjugate stabilizers.
A quaternion is a formal sum of a scalar and a 3D vector, aka real and imaginary parts. The scalar and vector components of the product of two vectors is $\mathbf{uv}=-\mathbf{u}\cdot\mathbf{v}+\mathbf{u}\times\mathbf{v}$, so to multiply two arbitrary quaternions you would FOIL and then use this formula. This entails a lot: the sqrts of $+1$ are $\pm1$, the sqrts of $-1$ are the unit vectors, two quaternions commute iff their vector parts are parallel, two vectors anticommute (meaning $\mathbf{vu}=-\mathbf{uv}$) iff they are perpendicular. The conjugate of a quaternion has the vector part negated, and the natural inner product is $\langle x,y\rangle=\mathrm{Re}(\overline{x}y)$, which basically means $\{1,\mathbf{i},\mathbf{j},\mathbf{k}\}$ is an orthonormal basis. Nonreal quaternions have unique polar forms $p=re^{\theta\mathbf{v}}$ where $r>0$ is the magnitude, $0<\theta<\pi$ is a convex angle, and $\mathbf{v}$ is a unit vector. You can expand this with $\cos$ and $\sin$ with Euler's formula. The effect of conjugation $p\mathbf{x}p^{-1}$ on 3D vectors $\mathbf{x}$ is to rotate it around $\mathbf{v}$ by an angle of $2\theta$.
As a group action this gives a double covering $S^3\to\mathrm{SO}(3)$ with kernel $S^0=\{\pm1\}$, which means $S^3$ is a spin group $\mathrm{Spin}(3)$. By using both left and right multiplication simultaneously, every pair $(p,q)$ of quaternions yields a 4D rotation $x\mapsto px\overline{q}$ (treating quaternions themselves as a 4D Euclidean space). This gives a double covering $S^3\times S^3\to\mathrm{SO}(4)$, so we can say that $S^3\times S^3\cong\mathrm{Spin}(4)$.
Here's where your formula for the fibers comes from: given two nonparallel unit vectors $\mathbf{u}$ and $\mathbf{v}$, we can classify the set of unit quaternions ("versors") $p$ for which $p\mathbf{u}p^{-1}=\mathbf{v}$. Its just a coset of the stabilizer $\mathrm{Stab}(\mathbf{u})$, which is the unit circle subgroup of $\mathbb{R}[\mathbf{u}]\cong\mathbb{C}$. To know which stabilizer, it suffices to exhibit a single coset representative. We could use $p=e^{\phi\mathbf{w}}$ where $\mathbf{w}$ is perpendicular to $\mathbf{u}$ and $\mathbf{v}$ and $\phi$ is half the angle between $\mathbf{u}$ and $\mathbf{v}$; with the half-angle formula this would yield
$$ p = \sqrt{\frac{1+\mathbf{u}\cdot\mathbf{v}}{2}} + \sqrt{\frac{1-\mathbf{u}\cdot\mathbf{v}}{2}}\frac{\mathbf{u}\times\mathbf{v}}{\|\mathbf{u}\times\mathbf{v}\|}, $$
which you could write in other formats if you wanted. Or, you could use an axis $\mathbf{z}$ halfway between $\mathbf{u}$ and $\mathbf{v}$ and a $180^\circ$ rotation, corresponding to
$$ p = \exp\left(\frac{\pi}{2}\mathbf{z}\right)=\mathbf{z}=\frac{\mathbf{x}+\mathbf{y}}{\|\mathbf{x}+\mathbf{y}\|}.$$
Any other choice of $p$ uses an axis on the great circle of $S^2$ containing $\mathbf{w}$ and $\mathbf{z}$. Note the axis $\pm\mathbf{w}$ is perpendicular to the axis $\pm\mathbf{z}$. Also note the left coset $p\mathrm{Stab}(\mathbf{u})$ is the right coset $\mathrm{Stab}(\mathbf{v})p$.
Proposition. Every great circle of $S^3$ is a coset of a circle subgroup.
[ Suppose $C$ is a great circle of $S^3$. Pick two perpendicular nonreal quaternions $p$ and $q$ from it. Then $\mathbf{u}=p^{-1}q$ is a unit vector, and $C$ is the coset $p\{\exp(\theta\mathbf{u})\}$. ]
Note you're talking about oriented 2D subspaces, hence oriented great circles. To accommodate for that in the above, replace circle subgroup with one-parameter circle subgroup.
It's a nice idea to use all Hopf fibrations (over all choices of basepoints) simultaneously. There's a pretty way to package this all together using the parallelizability of $S^3$. On the one hand, there's a kind of "currying" map $S^3\times S^2\to S^2\times S^2$ given by $(p,\mathbf{v})\mapsto (p\mathbf{v}p^{-1},\mathbf{v})$, where a Hopf map is applied to the first component and the second component determines which Hopf map is applied. On the other hand, $S^3\times S^2 = \mathrm{UT}S^3$ is the unit tangent bundle of $S^3$ via the correspondence $(p,p^{-1}q)\leftrightarrow (p,q)$, and each pair $(p,q)$ of perpendicular quaternions induces an oriented 2D subspaces. The projections $S^3\times S^2\to S^2\times S^2$ and $\mathrm{UT}S^3\to\widetilde{\mathbb{G}}_2\mathbb{R}^4$ both have circular fibers, which correspond to each other via $S^3\times S^2\simeq\mathrm{UT}S^3$.
In summary, we have a bundle isomorphism
$$ \begin{array}{ccccc}
S^1 & \longrightarrow & S^3\times S^2 & \longrightarrow & S^2\times S^2 \\
\updownarrow & & \updownarrow & & \updownarrow \\
S^1 & \longrightarrow & \mathrm{UT}S^3 & \longrightarrow & \widetilde{\mathbb{G}}_2\mathbb{R}^4
\end{array} $$
for which on the right half we can element-chase
$$ \begin{array}{ccc}
(p,p^{-1}q) & \mapsto & (qp^{-1},p^{-1}q) \\
\updownarrow & & \updownarrow \\
(p,q) & \mapsto & \mathrm{span}\{p,q\}
\end{array} $$
In other words, $\mathrm{span}\{p,q\}\in\widetilde{\mathbb{G}}_2\mathbb{R}^4$ (with orientation interpreted so $q$ is a positive right angle from $q$) corresponds to the pair of unit vectors $(qp^{-1},p^{-1}q)\in S^2\times S^2$. Note for any 2D subspace $\Pi$ there is an $S^1$s worth of choices of ordered orthonormal bases $\{p,q\}$, but $qp^{-1}$ and $p^{-1}q$ do not depend on this choice so this is all well-defined.
Here's another way to see $\widetilde{\mathbb{G}}_2\mathbb{R}^4\simeq S^2\times S^2$.
In general, we can embed $\widetilde{\mathbb{G}}_k\mathbb{R}^n\hookrightarrow\Lambda^k\mathbb{R}^n$ via $\mathrm{span}\{v_1,\cdots,v_k\}\mapsto v_1\wedge\cdots\wedge v_k$ (using ordered orthonormal bases, this is well-defined). The Hodge-star operator $\star:\Lambda^k\mathbb{R}^n\leftrightarrow\Lambda^{n-k}\mathbb{R}^n$ linearly extends orthogonal-complmentation between $\widetilde{\mathbb{G}}_k\mathbb{R}^n\leftrightarrow\widetilde{\mathbb{G}}_{n-k}\mathbb{R}^n$. In the event of $n=2k=4$, the star operator as $\pm1$-eigenspaces, each 3D. The space $\Lambda^2\mathbb{R}^4$ inherits an inner product from $\mathbb{R}^4$, in which case $\Lambda^2=\Lambda^2_+\oplus\Lambda^2_-$ is an orthogonal direct sum, with $a\wedge b\leftrightarrow a\wedge b\pm c\wedge d$ for every ordered orthonormal basis $\{a,b,c,d\}$. We can verify $\widetilde{\mathbb{G}}_2\mathbb{R}^4$ corresponds to $S^2\times S^2$ wrt this.
It is a standard fact that $\mathfrak{so}(n)\cong\Lambda^2\mathbb{R}^n$, where $a\wedge b$ represents the right-angle rotation from $a$ to $b$ (assuming $a,b$ are perpendicular unit vectors) in the $ab$-plane and acts as the $0$ map in the complement of the $ab$-plane. The $\star$-eigendecomposition basically corresponds to $\mathfrak{so}(4)\cong\mathfrak{so}(3)\oplus\mathfrak{so}(3)$, which is the infinitessimal version of $\mathrm{Spin}(4)\cong\mathrm{Spin}(3)\times\mathrm{Spin}(3)$.
In conclusion, $\widetilde{\mathbb{G}}_2\mathbb{R}^4\simeq S^2\times S^2$ is like the infinitessimal version of $\mathrm{Spin}(4)\cong S^3\times S^3$.
