If you just want to have a lower-bound on the time of classical existence, it is not too hard, though the answer will be not as sharp as the case of the Burgers' equation.

We will assume that $f_n$ is not just bounded, but also strictly positive; this is to ensure hyperbolicity of the system.

Then the function $\sqrt{f_n}$ is well-defined, and you can set $F_n(s) = \int_0^s \sqrt{f_n(\sigma)} ~\mathrm{d}\sigma$. Since $\sqrt{f_n}$ is strictly positive, we have that $F_n: \mathbb{R}\to \mathbb{R}$ is a diffeomorphism to its image, and is in particular invertible.

Your equation can be rewritten as $$ \partial_t p = \sqrt{f_n(r)} \partial_x F_n(r) , \qquad \frac{1}{\sqrt{f_n(r)}} \partial_t F_n(r) = \partial_x p $$ So if you redefine $\rho_n = F_n(r)$, and $g_n(\sigma) = \sqrt{f_n(F_n^{-1}(\sigma))}$, you have that the system becomes $$ \partial_t p = g_n(\rho_n) \partial_x \rho_n , \qquad \partial_t \rho_n = g_n(\rho_n) \partial_x p $$ Doing one final change of variables setting $u_n = \rho_n + p$ and $v_n = \rho_n - p$, you get that the equations satisfied are $$ \partial_t u_n = g_n(u_n - v_n) \partial_x u_n $$ and $$ \partial_t v_n = - g_n(u_n - v_n) \partial_x v_n $$

(Remark, the subscript $n$ here is to remind ourselves that these values, even at $t = 0$, depend on the choice of $f_n$.)

(Remark, finding $u_n$ and $v_n$ as above is basically doing the method of Riemann invariants by hand.)

Conclusion 1 Using the method of characteristics, within the classical life-span we have that $u_n$ and $v_n$ are bounded by their initial values. As all the variable transformations we made are invertible, you can use this to get a priori bounds on the $(p,r)$ variables in $L^\infty$ based on their initial values (no need to consider their derivatives).

Conclusion 2 Taking spatial derivatives, and using $D_t^\pm$ to denote the material derivatives $\partial_t \pm g_n(u_n - v_n) \partial_x$, we have $$ D_t^- u_n' = g'(u_n-v_n) [ (u_n')^2 - (u_n')(v_n') ] $$ and $$ D_t^+ v_n' = g'(u_n-v_n) [ (v_n')^2 - (u_n')(v_n') ] $$

Now, by our first conclusion we know that at any point where the classical solution is well-defined, $$ \min u_n(0,y) - \max v_n(0,y) \leq u_n(t,x) - v_n(t,x) \leq \max u_n(0,y) - \min v_n(0,y) $$ On this interval necessarily $|g'|$ is bounded, say by the number $\mathring{g} > 0$. Then the differential equations for $u_n'$ and $v_n'$ implies the following integral inequality: let $M(t) = \max_y |u_n'(t,y)| + \max_y |v_n'(t,y)|$ then

$$ M(t) \leq M(0) + \mathring{g} \int_0^t M(s)^2 ~\mathrm{d}s $$

A standard bootstrap argument shows that on $$ t \in [0, \frac{1}{C^2 \mathring{g} M(0)}] \tag{time}$$ we can guarantee that $M(t) \leq 2 M(0)$ and shock cannot form. Notice that the upper bound is essentially determined by the following pieces of data:

• The maximum of the derivative of $g_n$ on the interval allowed by the initial values of $u_n$ and $v_n$
• The maximum of the derivatives of $u_n$ and $v_n$

showing that this is qualitatively very similar to the Burgers' case you quoted.

Conclusion 3 To answer your second question: the answer is not quite. What we have are two-fold:

1. Up to the lower bound of the classical lifespan given by (time), indeed you can bound the $W^{k,\infty}$ norm of the solution by the $W^{k,\infty}$ norm of the initial data.
2. Let $T$ be any time such that the classical solution exists on time $[0,T]$, then you can bound the $W^{k,\infty}$ norm of the solution by the $W^{k,\infty}$ norm of the initial data together with the quantity $\max_{[0,T]} M(t)$ where $M(t)$ is defined as above in conclusion 2.

(In the case of Burgers, the two coincide because of our sharp estimate of $M(t)$; in the general cases, this is essentially the best we can do.)

Both claims follow by looking at higher derivatives of $u_n$ and $v_n$, which satisfy

$$ D_t^- u_n^{(k)} = U_k u_n^{(k)} + \hat{U}_k v_n^{(k)} + \tilde{U}_k $$

and

$$ D_t^+ v_n^{(k)} = V_k v_n^{(k)} + \hat{V}_k u_n^{(k)} + \tilde{V}_k $$

where $U_k, V_k$ and $\hat{U}_k, \hat{V}_k$ are expressions depending only on $g'(u_n - v_n)$, $u'_n$, and $v'_n$, while $\tilde{U}_k$ and $\tilde{V}_k$ are expressions depending on

• $g^{(j)}(u_n - v_n)$ for $1 \leq j \leq k$
• $u^{(j)}_n$ for $1 \leq j \leq k-1$
• $v^{(j)}_n$ for $1 \leq j \leq k-1$

and so inductively applying Gronwall's inequality will give you all the desired bounds.