Define a mean algebra to be a set $S$ with an binary operation $M$ satisfying (1), (2), and (4). We can define $M(a,b,c,d)=M(M(a,b),M(c,d))$ and this will depend only on the multiset $\{a,b,c,d\}$. More generally, we can think of $M$ as an operation defined on multisets of size $2^n$ for any $n>0$ (and this is well-defined by an easy induction on $n$ using (2) and (4)). The free mean algebra on two generators can be identified with the dyadic rationals between $0$ and $1$ with the arithmetic mean (the generators being $0$ and $1$); freeness of this algebra is easy to see when you write $k/2^n$ as $M(0,0,\dots,0,1,1,\dots,1)$ (with $k$ $0$s and $2^n-k$ $1$s). Denote this algebra by $Q$. Write $I$ for the algebra $[0,1]$ with the arithmetic mean; this can be thought of as a sort of completion of $Q$.
Suppose $A$ is a mean algebra with underlying set $[0,1]$ that further satisfies (5), (7), and (8). There is a unique mean-preserving map $f_0:Q\to A$ satisfying $f_0(0)=0$ and $f_0(1)=1$, and it follows easily from (5) and (7) that it is injective and order-preserving. The inverse of $f_0$ extends uniquely to an order-preserving surjection $g:A\to I$ (every element of $A$ defines a Dedekind cut in $Q$ via $f_0$). By (7) and (8) $g$ will also be mean-preserving, and so (5) implies $g$ is injective. We thus conclude that $A$ is isomorphic to $I$ as an ordered mean algebra, and furthermore this isomorphism is unique.
Now suppose $A$ is a mean algebra with underlying set $\mathbb{R}$ satisfying (5), (7), and (8). Every compact subinterval of $A$ these has a unique order-preserving isomorphism to $I$, and by uniqueness we can glue these isomorphisms together to get an order-preserving isomorphism between $A$ and an open subinterval of $\mathbb{R}$ with the arithmetic mean. Note that there are actually 4 distinct isomorphism classes of such structures (with respect to both the order and the mean): $(-\infty,\infty)$, $(-\infty, 0)$, $(0,\infty)$, and $(0,1)$. In particular, this gives an affirmative answer to Q4 (without assuming (6)). Furthermore, the isomorphism in question is clearly unique up to composition with affine maps $\mathbb{R}\to\mathbb{R}$.
As for your desire to only assume (3) and not (4), the only place where I used (4) is in asserting that $Q$ is free. Let $F$ be the free algebra on $\{0,1\}$ assuming (3) instead of (4); we wish to prove that the canonical map $F\to Q$ is an isomorphism. Every element of $F$ can be represented as a full binary tree of some height $n$ where each of the $2^n$ leaves is labelled by $0$ or $1$, and at each juncture of the tree we apply $M$. It suffices to show that if two such trees have the same number of leaves that are 1, then they represent equal elements of $F$. We prove this by induction on $n$; the cases $n\leq 1$ are trivial.
Let $x\in F$ be represented by such a tree of height $n>1$ that has $i$ leaves that are $1$. WLOG $i\leq 2^{n-1}$ (otherwise we can just swap the roles of $0$ and $1$). Let $y$ be represented by a tree that looks the same as $x$'s, except that all the $0$s are on the left and all the $1$s are on the right. For instance, if $n=i=3$ and $$x=M(M(M(1,0),M(0,1)),M(M(0,1),M(0,0))),$$ then $$y=M(M(M(0,0),M(0,0)),M(M(0,1),M(1,1))).$$ It suffices to prove that $x=y$, since $y$ depends only on $n$ and $i$. Furthermore, by induction, it suffices to prove that $y=M(z,w)$ where $z$ and $w$ each have as many $1$s as the left and right halves of $x$, respectively (since $z$ and $w$ can then be transformed to look the same as the two halves of $x$). Let $j$ be the number of $1$s in the left half of $x$; WLOG $j\leq 2^{n-2}$ (if not, switch the two halves of $x$). By the induction hypothesis, we can write the right half of $y$ as $M(b,c)$, where $b$ has $j$ $1s$. Also, the left half of $y$ is $M(a,a)$, where $a$ is the tree of height $n-2$ consisting entirely of $0$s. We thus have $y=M(a,M(b,c))$, and so by (3) we can rewrite it as $y=M(M(a,b),M(a,c))$. Setting $z=M(a,b)$ and $w=M(a,c)$, the proof is complete.
As a final remark, I'm not sure what happens when you additionally assume (6). Certainly most functions $f$ as in your Q4 will not give rise to a mean satisfying (6); the only ones I know of that do are $f(x)=x^p$ for $p\neq0$ and $f(x)=\log x$, up to composition with affine maps. Note that $f(x)=\log x$ (corresponding to the geometric mean) can be thought of as the $p=0$ case; indeed, the mean obtained from it is equal to the limit of the $x^p$ means as $p\to 0$, and these functions $f$ (including their compositions with affine maps) are exactly the solutions of the differential equation $(xf''/f')'=0$. Perhaps if you assume $f$ is sufficiently differentiable you could prove it must be of this form by differentiating the functional equation you get from associativity.
EDIT: Here's a proof that if (6) holds, then $f(x)$ must be of the form $ax^p+b$ or $a\log x+b$ and hence $M$ is either $((x^p+y^p)/2)^{1/p}$ or $\sqrt{xy}$. Suppose that (6) holds and $f:\mathbb{R}_+\to U\subseteq\mathbb{R}$ is an isomorphism from $M$ to the arithmetic mean on some open interval $U$. For $c>0$, consider the map $x\mapsto f(cf^{-1}(x))$. This is a mean-preserving automorphism of $U$, so it must be of the form $x\mapsto A(c)x+B(c)$ for some constants $A(c)$ and $B(c)$ depending on $c$. That is, we have $f(cx)=A(c)f(x)+B(c)$.
Write $r=A(e)$; then it is easy to show that if $c=e^q$ for $q\in\mathbb{Q}$, we have $A(c)=r^q$. By continuity, it follows that $A(c)=r^{\log c}=c^p$ for all $c$, where $p=\log r$. Suppose that $p=0$, so $f(cx)=f(x)+B(c)$. Then $B(cd)=B(c)+B(d)$ and it follows by continuity that $B(c)=a\log c$ for some $a$. Setting $b=f(1)$, we then get $f(c)=a\log c +b$, as desired.
Now suppose that $p\neq 0$. Let $s=B(e)$; then by induction we have $$B(e^n)=s\frac{r^{n+1}-1}{r-1}$$ for $n\in \mathbb{N}$. But by the same argument using $e^{1/m}$ instead of $e$ for some integer $m>0$, we must have $$B(e^n)=B(e^{1/m})\frac{r^{n+1}-1}{r^{1/m}-1}.$$ It follows that the first formula is also valid when $n=1/m$, and it is now easy to show that it is valid for all rational $n$. By continuity, it holds for all real $n$, and we can rewrite it to get $$B(c)=s\frac{e^pc^p-1}{r-1}=tc^p+b$$ for some constants $t$ and $b$. Thus we have $f(cx)=c^pf(x)+tc^p+b$, and setting $x=1$ and $a=f(1)+t$ we get $f(c)=ac^p+b$, as desired.