For the second part, if there exists a simultaneous solution there must exists constants $\lambda,\mu$ such that
$ r = \lambda a + \frac{a\wedge b}{|a|^2} = \mu c + \frac{c\wedge d}{|c|^2} $
by the first part of your question. Then
$ r\wedge c = d \implies \lambda a \wedge c + \frac{1}{|a|^2}(a\wedge b) \wedge c = d \tag{1}$
and
$ r\wedge a = b \implies -\mu a \wedge c + \frac{1}{|c|^2}(c\wedge d)\wedge a = b $
We consider (1). The second equation can be solved analogously.
Rewrite (1) as
$ \lambda a\wedge c = d - \frac{1}{|a|^2}(a\wedge b)\wedge c $
we see that to have a solution $\lambda$ we need that the vector $a\wedge c$ and the vector $d - \frac{1}{|a|^2}(a\wedge b)\wedge c$ are collinear, in which case we can "divide by $a\wedge c$" (since we assumed that $a\wedge c \neq 0$) and get the correct factor $\lambda$. Now, we know that $a\wedge c$ is orthogonal to both $a$ and $c$. Since we are working in three dimensions, to show that the right hand side is collinear with $a\wedge c$ it suffices to show that it is also orthogonal to both $a$ and $c$.
Now recall the vector triple product formula which gives $(a\wedge b)\wedge c = -(c\cdot b) a + (c\cdot a)b$.
Taking the dot products we require
$ 0 = d\cdot c - \frac{1}{|a|^2} \left[ (c\cdot a)(b\cdot c) - (c\cdot b)(a\cdot c)\right] \tag{2}$
and
$ 0 = d\cdot a - \frac{1}{|a|^2} \left[ (c\cdot a)(b\cdot a) - (c\cdot b)(a\cdot a)\right] \tag{3}$
Now, under the assumption that $c\cdot d = 0$, equation (2) holds automatically. Under the assumption that $b\cdot a = 0$, equation (3) holds if and only if $ 0 = a\cdot d + b\cdot c$ which is the condition you desire.