Consider the mixed formulation of the Poisson/Darcy system for a region $\Omega$:
$\alpha \mathbf{v} + \nabla p = f \\ \mathrm{div}[\mathbf{v}] = 0 $
with the boundary conditions
$\mathbf{v}\cdot \mathbf{n} = v_n$ on $\Gamma_v$,
$p = p_0$ on $\Gamma_p$
For the weak form: $a(\mathbf{v}, \mathbf{w}) := \int_\Omega \mathbf{v} \cdot \mathbf{w}$,
while $b_1(\mathbf{w},p) := \int_\Omega p\;\mathrm{div}[\mathbf{w}]$ or
$b_2(\mathbf{w},p) := \int_\Omega \nabla p\cdot\mathbf{w}$.
Depending on the choice for $b_1$ or $b_2$, the solution spaces can be chosen as:
1) $\mathbf{v} \in H(div; \Omega)^n $ while $p \in L^2(\Omega)$ for the choice $b_1$. One can choose for this the BDM elements for the velocity and piecewise constant pressures for each element and get a stable solution. In this formulation, the condition on velocity is essential while the boundary condition for pressure is natural.
2) $\mathbf{v} \in L^2(\Omega)^n, p \in H^1_0(\Omega)$ for the choice $b_2$. In this case, the velocity can be approximated by the discontinous Raviart-Thomas elements of polynomial order $k$ and the pressure by Lagrangian elements of polynomial order $k+1$, to get a stable solution. In this formulation, the condition on velocity is natural while the boundary condition for pressure is essential.
Q1. Is one of 1), 2) to be preferred over the other choice? Will the solutions obtained differ ?
Q2. However, $H^1_{0}(\Omega)^n \subset H(div;\Omega)^n$. So suppose I choose $\mathbf{v} \in H^1_0(\Omega)^n, p \in L^2(\Omega)$ for the choice $b_1$ and use the standard Taylor-Hood element, enforcing velocity boundary condition as essential and the pressure boundary condition as natural. Is this a stable approximation ?
PS: I am unable to create the tags "Sobolev-spaces", "mixed-formulations". If you agree that they are apt for this question, could you please create them.