The same content is exposed in a section of Ref. (1) on the Buckley-Leverett equation as follows (p. 352-353):
We only need to solve for $q$, the water saturation, and so we can solve a scalar conservation law
$$ q_t + f(q)_x = 0 $$
with the flux $f(q)$ given by [
$$ f(q) = \frac{q^2}{q^2 + a (1-q)^2} . \tag{16.2} $$
Here $a<1$ is a constant]. Note that this flux is nonconvex, with a single inflection point.
and
Now consider the Riemann problem with initial states $q_l = 1$ and $q_r = 0$. By following characteristics, we can construct the triple-valued solution shown in [Fig. 4.7(a) of OP]. Note that the characteristic velocities are $f'(q)$, so that the profile of this bulge, seen here at time $t$, is simply the graph of $t f'(q)$ turned sideways.
Thus, we need to solve a Riemann problem in the nonconvex case. We consider a self-similar solution of the form $q(x,t) = v(\xi)$ with $\xi = x/t$. Using the chain rule, we have $q_t = v'(\xi)\xi_t$ and $q_x = v'(\xi)\xi_x$ with $\xi_t = -\xi/t$ and $\xi_x = 1/t$. Injecting the previous Ansatz in the quasi-linear PDE $q_t + f'(q)q_x = 0$ satisfied by the water saturation leads to
$$
\left(f'(v(\xi)) - \xi\right) {v'(\xi)} = 0\, ,
$$
which nontrivial solutions satisfy $f'(q(x,t)) = x/t$. Hence, the bulge at some time $t$ is obtained by plotting $x = t f'(q)$ for $q$ in $[0,1]$. This is illustrated below, where $a=1/2$ and $t=1$ (abscissas $q$, ordinates $x$):
It remains to rotate the figure to obtain $q$ as a function of $x$. This topic is also well-explained in an interactive book, more precisely §1.6 and §5.1 on nonconvex conservation laws.
The equal area rule results from conservation, more precisely from the fact that a discontinuous weak solution must have the same integral as the smooth multivalued solution deduced from the characteristics. This is exactly the way it is applied in Fig. 4.7(b) of OP.
(1) R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002, doi:10.1017/CBO9780511791253
Best Answer
Let us consider the scalar one-dimensional case $u_t + f(u)_x = 0$ for sake of simplicity, as described in Section 3.8 of (1).
An entropy function $\eta$ is a convex function of $u$ which satisfies the conservation law $\eta(u)_t + \psi(u)_x = 0$, where $\psi$ is the entropy flux. In the definition, neither the flux $f$ nor the entropy flux $\psi$ is convex, but only $\eta$ is convex. This is a conventional choice with important consequences, "though one may consider a more general framework" (2). If $\eta$ is chosen convex, then the vanishing viscosity weak solution -- i.e., the physically relevant weak solution -- makes the entropy decrease in the weak sense. This leads to the definition of the entropy solution, which is a weak solution $u$ such that $\eta(u)_t + \psi(u)_x \leq 0$ weakly for all convex entropy function $\eta$.
As such, the definition of an entropy function $\eta$ is not related to the thermodynamical notion of entropy. Firstly, it decreases over time, while the entropy increases in thermodynamics. However, there may be a link with the thermodynamical entropy for some particular systems. In the case of Eulerian gas dynamics, a natural mathematical entropy function is defined with respect to the thermodynamical entropy (cf. Section II.2 of (2)).
(1) R.J. LeVeque, Numerical Methods for Conservation Laws. Birkhäuser, 1992. doi:10.1007/978-3-0348-8629-1
(2) E. Godlewski, P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws. Springer, 1996. doi:10.1007/978-1-4612-0713-9