Density and area in 1D conservation of momentum and nozzle design I'm following this derivation from the NASA website explaining the design of converging/diverging nozzles in rocket design, specifically the derivation in image form at the top of the page:
https://www.grc.nasa.gov/www/k-12/airplane/nozzled.html
In the end we arrive at an identity which states that when the Mach number is below 1, increase in area causes velocity to decrease, and the Mach number is over 1, increase in area cause flow velocity to increase.
One of the identities used in the derivation is the 1D conservation of momentum, which is derived here:
https://www.grc.nasa.gov/www/BGH/conmo.html
The derivation of the conservation of momentum assumes that the flow is incompressible, and therefore has a constant density. But how can it be therefore used to derive the area/velocity relationship (first link), which clearly needs to hold for both compressible and incompressible flows? It is even stated in the article:

For subsonic (incompressible) flows, the density remains fairly constant, so the increase in area produces only a change in velocity. But in supersonic flows, there are two changes; the velocity and the density.

I understand that this derivation is very simplified; but using an identity that seems to hold only for incompressible flows for a derivation dealing with both incompressible and compressible flows seems wrong.
EDIT: An additional question came into my mind: In the conservation of momentum derivation, why is it reasonable to also assume that the area is constant as well? This also seems to contradict everything in the nozzle design page as well as the fact that the pressure and velocity are changing (how can we assume that pressure and velocity change but the area does not?)
 A: The Conservation of Momentum page does not assume the incompressibility and, consequently, the equation derived holds for supersonic flows.
To clarify: up to and including the equation
$$
- \frac{\Delta p }{\Delta x} \cdot \Delta x \cdot A = m \left( \frac{\Delta u} {\Delta x}\right) \cdot \frac{\Delta x }{ \Delta t},\tag{1}
$$
density does not come up at all, we are dealing with the change of momentum of the small fluid volume as it flows along its single direction. Density inside this fluid element could, of course, change as it moves but its mass does not, and so only velocity derivative  could appear in the r.h.s. of (1).
Density of course, will appear at a next step, when we replace $m$ with $\rho \Delta x \,A$, but since we are keeping only the leading order of contributions, whether it is constant or not here does not matter.
After simplifying the result the remaining step is to go from difference relations to derivatives:
$$
-\frac{dp}{dx}=\rho u \frac{du}{dx}.
$$
This is in perfect agreement with Euler equation (taken from here)
$$
\frac{\partial \mathbf{ u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + \frac {\nabla p}{\rho} = \mathbf{g},
$$
if we assume stationary ($\frac\partial{\partial t}=0$) and one-dimensional ($\nabla=\frac d{dx}$) flow.
Your second question (why the area is assumed constant?) has similar answer: it isn't. The equation would remain valid as long as the flow remains oriented in one direction. 
Area first appears in the l.h.s.:
$$
F=- [(p \, A)_2 - (p \, A)_1] =-\left[\left(p + \frac{\Delta p }{\Delta x} \Delta x\right)  A - p \, A\right].
$$
Even if we assume that the surface surrounding the fluid has different $A_1$ and $A_2$ the result would remain the same in the first order of $\Delta x$.
The reason for this is simple: if we can consider pressure gradient $\nabla p$ approximately constant in the region of our small fluid volume then the force
 could be written as 
$$F=-\frac{dp}{dx} \cdot \mathrm{Vol}(m),$$
where Vol(m) is the volume of  (small) fluid element. And the geometry of that fluid element does not really matter, as long as it remains physically small. 
This is just like Archimedes' principle in hydrostatic (in fact the latter is simply a particular form of pressure gradient $dp/dx=\rho g$). And derivation of Archimedes' principle (here's a link on Physics SE) is also applicable in our case: we can approximate the volume of a fluid element with large amount of infinitesimally small 'straws' along the direction of pressure gradient and each straw would provide contribution $-(dp/dx) l\,\delta A$, where $l$ is the length of a given straw and $\delta A$ is area of its crossection.
So if we have the flow which is slightly expanding so that the area $A$ is a (slowly changing function of $x$) then the force on our fluid element could diverge from $(d p/d x) \Delta x \, A $ in three ways: (1) value of pressure gradient changes along the element, (2) direction of it changes along the element, (3) $A$ is not constant in an element. In each of the cases the error in the expression for the force is proportional to $(\Delta x)^2$ and so would disappear once we make the element sufficiently small. For example for (3):
$$
|\delta F_3| \le |A_2-A_1||p_2-p_1|\approx \left|\frac{dA}{dx}\right|\Delta x \left|\frac{dp}{dx}\right| \Delta x \sim (\Delta x)^2,
$$
A: 
But how can it be therefore used to derive the area/velocity relationship (first link), which clearly needs to hold for both compressible and incompressible flows?

Momentum conservation holds regardless of the compressibility of a fluid.  I think you are confusing a specific form of the continuity equation (i.e., often one uses a mass flux continuity equation) with a conservation of momentum.

An additional question came into my mind: In the conservation of momentum derivation, why is it reasonable to also assume that the area is constant as well?

I do not think the area is assumed to be constant (if it is assumed, then you are correctly confused because this is wrong).  You can follow the simple derivations, starting from an energy flux conservation equation, of how the acceleration of the flow depends upon the sign of dA by looking at the Wikipedia page for isentropic nozzle flow.  The end result is an expression of the fluid speed, $u$, versus the nozzle area, $A$, given by:
$$
\frac{du}{u} \ \left( M^{2} - 1 \right) = \frac{dA}{A} \tag{0}
$$
where $M$ = Mach number.  From Equation 0, one can see the following list:


*

*for $M$ < 1 and $dA$ > 0 $\rightarrow$ $du$ < 0 (i.e., deceleration);

*for $M$ < 1 and $dA$ < 0 $\rightarrow$ $du$ > 0 (i.e., acceleration);

*for $M$ > 1 and $dA$ > 0 $\rightarrow$ $du$ > 0; and

*for $M$ > 1 and $dA$ < 0 $\rightarrow$ $du$ < 0.



This also seems to contradict everything in the nozzle design page as well as the fact that the pressure and velocity are changing (how can we assume that pressure and velocity change but the area does not?)

Again, I think you misunderstood something.  I do not see a statement that the conservation of momentum requires $dA$ = 0, which it does not.
More detailed approach
Much of the following is explained in more detail in the answer found at:
https://physics.stackexchange.com/a/137842/59023
and is largely taken from Whitham [1999].
The Mach number can be written as a function of the area of the channel given by:
$$
\frac{ 1 }{ A } \frac{ d A }{ d M } = -g\left( M \right) \tag{1}
$$
where g(M) is given by:
$$
\begin{align}
g\left( M \right) & = \frac{ M }{ M^{2} - 1 } \left( 1 + \frac{ 2 }{ \gamma + 1 } \frac{ 1 - \mu^{2} }{ \mu }  \right) \left( 1 + 2\mu + \frac{ 1 }{ M^{2} } \right) \tag{2a} \\
\mu^{2} & = \frac{ \left( \gamma - 1 \right) M^{2} + 2 }{ 2\gamma M^{2} - \left( \gamma - 1 \right) } \tag{2b}
\end{align}
$$
where $\gamma$ = ratio of specific heats, $M$ = Mach number.  We can then rewrite Equation 1 (using the chain rule) to find:
$$
g\left( M \right) \frac{ d M }{ d x } + \frac{ 1 }{ A } \frac{ d A }{ d x } = 0 \tag{3}
$$
Note that $g\left( M \right)$ < 0 is satisfied on the interval $\left\{\sqrt{4 - \sqrt{13}}, 1 \right\}$.  Below $M$ ~ 0.628 the magnitude of $g\left( M \right)$ diverges toward $+ \infty$ near $M$ ~ 0.447.  Equation 3 is really only valid for values of $M$ near unity or larger than unity (probably breaks down for $M$ > 5 or 6 too), so when I say $M$ < 1 below, I am really implying 0.65 < $M$ < 1.
For $dA/dx$ > 0 one needs $g\left( M \right) \ dM/dx$ < 0, which means either the Mach number decreases as the fluid moves toward the nozzle exit (i.e., $dM/dx$ < 0) or the flow is subsonic, $M$ < 1 (i.e., $g\left( M \right)$ < 0).  If the sound speed were roughly constant (but definitely not increasing) from the inlet to the exit of the nozzle, then we would not expect $dM/dx$ < 0 unless the speed decreased.  Therefore, one would need $M$ < 1 to satisfy Equation 3 for $dA/dx$ > 0.
For $M$ > 1, then $g\left( M \right)$ > 0 for all $M$ which means that the sign of $dA/dx$ must be opposite of $dM/dx$ to satisfy Equation 3.
References
Whitham, G. B. (1999), Linear and Nonlinear Waves, New York, NY: John Wiley & Sons, Inc.; ISBN:0-471-35942-4.
