Is standard error propagation appropriate in this situation, and if so how, and why? I have a real-world problem which has left me puzzled in terms of how to approach it with any degree of rigour.  I'm not a statistician,, and although I can understand it enough to produce a vague answer, I'd like to actually understand the stats involved. 
In order to illustrate what I mean, I'll give a simplified example which outlines the problem.
Let's say I'm making a solution of salt. I use rock salt that is quoted as 85-95% sodium chloride.  In my solution I will add water, and let's say that the process gives us a range of 8%-12% rock salt in the solution.
My first reflex is to say that my end solution is therefore going to be 6.8%-11.4% sodium chloride (the product of the lower and upper bounds).  But this is such a wide range, and if you imagine the effects of this compounded (if I was to make solutions of solutions etc.) you would end up with ridiculously impractical end ranges that have no baring on reality.
So, it seems like the best thing to do is use "standard" error propagation.  Hopefully my maths is right here:
using:${\left( \frac {\sigma _x}{\bar x }\right)}^2= {\left(\frac {\sigma _y}{\bar y }\right)}^2+{\left(\frac {\sigma _z}{\bar z }\right)}^2$, I suppose I would say what I have is:
$(0.9±0.05) \cdot (0.1±0.02) = 0.09 \pm 0.09\cdot\sqrt{ {\frac{0.05}{0.9}}^2 + {\frac{0.02}{0.1}}^2}=0.09 \pm 0.0187$  i.e. giving me an end range of $0.071 < x < 0.109.$ (or 7.1%-10.9%)
What (I think) I've done here is assume that the real sample variation is normally distributed and that the upper and lower bounds represent 1 standard deviation from the mean of that distribution? 
But if it were the case instead that my errors instead represented say 2 sigma, or 3 sigma, would the above calculation still be valid? For example, if my sample ranges represent 2 standard deviations (i.e. 95% of all samples will fall within that range), and the same for the upper and lower bounds of the amount of water added, is it mathematically valid to do exactly the same calculation above, and just interpret that the end error (±0.019) therefore represents 2 standard deviations in the end sample?  If so, why is this valid, and if not, how should my calculation change? And if I assume that those initial inequalities represent, say, 1 sigma, how do I give an error for the final product that that represents 2 sigma, 3 sigma etc. etc. 
Bonus question: This one is a bit more complicated, so if you only answer the first part I'll still be really happy. Let's say I'm using a material that is definitely not normally distributed.  To give the extreme case, lets say instead that there is a source material which exists in nature in some (unknown) distribution of purity, but that this material is then divided up into grades (60-70%, 70-80%, 80-85%,85-90% etc.), based on a test with an accuracy much narrower than the ranges themselves.  So for the sake of simplicity, let's say the test is accurate enough that the error is neglible, and so functionally the distribution of whichever purity of material I use is not gaussian, but has a binary end. We don't know the distribution of the source material, so we essentially take the "graded" material as being just being a "flat" distribution with sharp edges (i.e. ___|‾|___ ).   Let's say, for the sake of argument, that my process can be best modelled this way too.  (Just as a thought experiment, let's say I add a random amount of water to each batch, and then for each batch I take a sample, evaporate the water, measure the solids content, and assign the batch to a subcategory of concentration ranges based on my measurement.)
So I have 2 flatly distributed ranges, and I'm taking the product of them.  What does the distribution of my end material look like, and how best should I express this as a value with an "error".  If I just completely ignored the flatness of the source data, and modelled it instead as 2 gaussian distributions, where the edges of the ranges are 2 sigma or 3 sigma, would I still end up with a "reasonable" end result for error bounds?
 A: 
What (I think) I've done here is assume that the real sample variation is normally distributed and that the upper and lower bounds represent 1 standard deviation from the mean of that distribution?

Not really. You actually didn't need to make any assumption on the underlying distribution nor did you. 

is it mathematically valid to do exactly the same calculation above, and just interpret that the end error (±0.019) therefore represents 2 standard deviations in the end sample?

No. Let's generalise this question a bit here and let's call the percentage of NaCL in the rock salt the RV $Y$ and the percentage of rock salt used in the mixture the RV $Z$. Now we essentially wish to understand the distribution of the product $X = YZ$. More specifically we really only care about the variance of this RV (or its sd) so we can understand how it changes when changing the variances (or sd's) of $Y$ and $Z$, which makes life easier as the pdf is a bit messier to derive.
The variance of the product of RV's is easy to solve for and is distribution independent. It is essentially what you are adding to the expectation of the product $E[YZ]$. As a proof (under independence assumptions):
$$\begin{align}
\operatorname{var}(X_1\cdots X_n) 
&= E[(X_1\cdots X_n)^2]-\left(E[X_1\cdots X_n]\right)^2\\
&= E[X_1^2\cdots X_n^2]-\left(E[(X_1]\cdots E[X_n]\right)^2\\
&= E[X_1^2]\cdots E[X_n^2] - (E[X_1])^2\cdots (E[X_n])^2\\
&= \prod_{i=1}^n \left(\operatorname{var}(X_i)+(E[X_i])^2\right)
- \prod_{i=1}^n \left(E[X_i]\right)^2
\end{align}$$
Note that I did not require any distribution assumptions (I guess other than defined first and second moments and independence between variables which I stated above). For two variables this reduces to:
$$ {\rm Var}(XY) = E(X^2Y^2) − (E(XY))^2={\rm Var}(X){\rm Var}(Y)+{\rm Var}(X)(E(Y))^2+{\rm Var}(Y)(E(X))^2$$
Which is essentially the square of what you are adding since you define the bounds as the plus minus of 1 sd.
Now knowing the variance of the product as a function of the variances of the marginals we can see how it changes with constant scaling. If we wish to determine with a scale of 2, the right hand side of our function becomes.
$$
\begin{align}
{2*\rm Var}(X)*2*{\rm Var}(Y)+{2*\rm Var}(X)(E(Y))^2+2*{\rm Var}(Y)(E(X))^2
&\neq 2*\rm Var(XY)
\end{align}
$$
Bonus Question: The "flat" distribution you're describing is commonly known as the uniform distribution, and like I said previously it wouldn't change anything as the variance of the product of random variables is distribution invariant.
Using an example for clarity: Let's say the percentage of NaCL is some random amount $Y$ which has error bounds ± $a$ where $a$ represents 1 standard deviation of $Y$. Likewise for the percentage of rock salt, lets call it $Z$ and it can vary with error bounds ± $b$ where $b$ represents 1 standard deviation of $Z$. Since they're both essentially concentration percentages, when we mix the two the final concentration of the mixture will essentially have a concentration amount of $Y \times Z$ which we can call $X$. Now obviously $X$ is a random variable, and one standard deviation of it's distribution is what you wrote before and what I proved above but for this examples sake it is essentially: 
$$sd(X) = \sqrt{a^2 b^2 + a^2 E(Z)^2 + b^2 E(Y)^2}$$
Now if we let the error bounds $a$ and $b$ now represent 2 standard deviations the question you asked was can we just use the formula above to represent two standard deviations in the product $X$. Let's see. The variance of the product is still given by the exact same formula, but since now $a$ and $b$ represent 2 sd's we must use $\frac{a}{2}$ and $\frac{b}{2}$ as 1 sd. So we get:
$$\begin{align}
sd(X) &= \sqrt{(\frac{a}{2})^2 (\frac{b}{2})^2 + (a/2)^2 E(Z)^2 + (b/2)^2 \times E(Y)^2}\\
&= \frac{1}{2}\sqrt{\frac{1}{2}a^2 b^2 + a^2 \times E(Z)^2 + b^2 \times E(Y)^2}
\end{align}$$
And if we double this to see what 2 standard deviations of the product look like we can see that it is not simply 2 times the orignial formula above, we will get:
$$
\begin{align}
 2 \times sd(X) &= \sqrt{\frac{1}{2}a^2 b^2 + a^2 \times E(Z)^2 + b^2 \times E(Y)^2}\\
&\neq \sqrt{a^2 b^2 + a^2 E(Z)^2 + b^2 E(Y)^2}
\end{align}
 $$
 So simply using the formula you provided does not scale when $a$ and $b$ represent values other than 1sd. There is an extra $\frac{1}{2}$ in front of the $a^2 b^2$ factor in this case, and more generally for when $a$ and $b$ represent $m$ standard deviations, there will be a factor of $\frac{1}{m}$.

and if not, how should my calculation change?

For $a$ and $b$ representing $m$ standard deviations, you will need to scale the calculation by a factor of:
$ 
\frac{\sqrt{\frac{1}{m}a^2 b^2 + a^2 \times E(Z)^2 + b^2 \times E(Y)^2}}{\sqrt{a^2 b^2 + a^2 E(Z)^2 + b^2 E(Y)^2}}
$
