Measuring $w$ is actually what I do for a living. The current best measurements put $w$ at $-1$ but with an uncertainty of $5\%$, so there's a little room for $w \ne -1$ models, but it's not big and getting smaller all the time. Indeed, we'd all be thrilled if, as measurements got more precise, $w \ne -1$ turns out to be the case because the $\Lambda$CDM model (where Dark Energy is just a scalar correction to General Relativity) is just horribly arbitrary and we'd much rather get evidence for something new and exciting.
You ask if models can reproduce a $w \approx -1$ and then ask about free parameters in such models. You're on the right track with asking about free parameters. There are some ways in which you can explain that $\Omega_\Lambda$ and $\Omega_m$ are fairly close to each other today (quintessence models do this with tracker models, inhomogeneity inducing corrections to the Friedmann equations can do this because the inhomogeneities get stronger as structure grows) while the $w \approx -1$ value is not a particularly unnatural value for many models (it is for some, and those are already ruled out). So you're in an odd situation where you could certainly imagine, say, a quintessence model with say $w \approx -1$, but you're not too tempted to because $\Lambda$CDM is just much simpler. To get out of this impasse, you're on the right track when you ask about free parameters (or loose ends as you call them). Quintessence for example would have a total of $4$, which represents the number of free parameters in your quintessence field that have an impact on the rate of evolution, whereas $\Lambda$CDM has only one. This is why, for the time being, we go with $\Lambda$CDM : our current probes can only constrain a single parameter, namely the current effective equation of state, and it is consistent with $\Lambda$CDM, which then plays the role of a sort of "placeholder" model given that it requires the least number of degrees of freedom. Current constraints on the evolution of Dark Energy properties are a joke, and you might as well consider them non existent. The image below is a confidence contour in the $w$ VS $w_a$ plane given the latest CMB, BAO, and supernova data. $w_a$ is a somewhat arbitrarily chosen parameter ($w(a) = w_0 + (1 - a) w_a$ where $a$ is the scale factor) but it gives you an idea of your probe's sensitivity to the evolution of $w$.
Future probes may be able to meaningfully constrain $w_a$, and if it turns out to be different than $0$, then that'll rule out $\Lambda$CDM and open up a whole new ball game. Remember, $\Lambda$CDM corresponds to a static $w_a = 0$. It'll be hard, since we'll need to go to much higher redshifts while maintaining similar levels of accuracy in our measurements, if not improvements because given its current properties (and whatever its evolution) Dark Energy is weaker in the past ($\Omega \propto a^{-3(1+w)}$, so while a $w \approx -1$ Dark Energy maintains a similar density, matter density, which corresponds to $w = 0$, increases as you look back in time). Being able to better constrain this evolution is the name of the game if you're interested in ruling out the current standard model.
The classification of possible explanations that your review suggests is indeed the best way to broadly categorize the proposed solutions.
A popular model which modifies the right hand side is quintessence. Quintessence employs a common solution (or common attempt at a solution) to many problems : When in doubt, throw in a scalar field. No known field satisfied the required properties, so it would most certainly lie outside the Standard Model of particle physics. In general, one can consider the impact on expansion of the addition of any fundamental field to the Lagrangian of particle physics. The only condition this field would have to satisfy for this field to be a candidate for dark energy is $w < -1/3$. For a quintessence field $\phi$ with a potential field $V(\phi)$, $p$ and $\rho$ are such that $w$ is given by :
\begin{eqnarray}
w_\phi &=& \frac{p}{\rho} & \mbox{ (by definition)} \\
&=& \frac{\frac{1}{2} \dot{\phi}^2 - V(\phi) - \frac{1}{6a^2} (\nabla \phi)^2}{\frac{1}{2} \dot{\phi}^2 + V(\phi) + \frac{1}{2a^2} (\nabla \phi)^2} & \mbox{ ($p$ and $\rho$ for quintessence fields)} \\
&=& \frac{\frac{1}{2} \dot{\phi}^2 - V(\phi)}{\frac{1}{2} \dot{\phi}^2 + V(\phi)} & \mbox{ (use spatial homogeneity to get rid of $\nabla \phi$ terms)} \\
&\stackrel{\text{$\dot{\phi} \rightarrow 0$}}{\approx}& -1 & \mbox{ (for a slow rolling field)}
\end{eqnarray}
The couple of $\nabla \phi$ terms in there disappear when you invoke a homogeneous universe (using the Comological Principle). Hence, many models with a $V(\phi)$ that yields a $w_\phi < -1/3$ have been considered simply by constructing the proper $V(\phi)$, and it is not hard nor particularly unnatural to construct models where $w \approx -1$ today, simply by having a slow rolling field where $\dot{\phi}$ is small and therefore $w_\phi \approx -1$. The difficulty lies in making such models distinguishable from the standard case of a $w = -1$ that is static in time, given our experimental insensitivity to the time evolution of the field. The field hasn't necessarily been in slow roll for the history of the universe. On the whole this is actually a whole class of models, and many models and couplings with matter can be invoked to explain all sorts of things about Dark Energy (like the rough similarity of $\Omega_\phi$ and $\Omega_m$ today, and why this happens when the field is slowly rolling etc etc). Again though, we're not yet sensitive to any of the observational implications of any of this. In the near future, we may be able to constrain 2 parameters of freedom with regards to the evolution of expansion, and not the 4 that are required to fully characterize a quintessence model, meaning that even then we will only be able to, at best, favor quintessence models over others, but it'll still be a stretch before we can definitively prove them.
Another right hand side explanation is k-essence, but this is essentially a modified version of quintessence where the kinetic energy term is not the canonical one, and rather than play around with $V(\phi)$ it is thought that it is the properties of this kinetic term that could lead to the required properties. Recall that our own kinetic term was $\frac{1}{2} \dot{\phi}^2 - \frac{1}{6a^2} (\nabla \phi)^2$, and that we got rid of it entirely by assuming $\nabla \phi$ homogeneity and that $\dot{\phi}$ was small. The idea for k-essence is not to get rid of these terms, but to come up with a new form for the kinetic energy such that the k-essence field will have to desired properties. The novelty here (relative to quintessence) is that we don't need a very low upper bound on the mass of field to make it slow rolling. But observationaly, this remains just as hard to prove because the properties of the field (your review gives a very comprehensive overview of this) leave open a number of degrees of freedom as well.
$f(R)$ models modify the left hand side. Here, again, it's best to think of this as a class of models. As the name indicates, $f(R)$ gravity works by generalizing the Einstein-Hilbert action :
$$
S[g] = \int \frac{1}{16 \pi G c^{-4}} R \sqrt{-|g_{\mu \nu}|} d^4x \rightarrow S[g] = \int \frac{1}{16 \pi G c^{-4}} f(R) \sqrt{-|g_{\mu \nu}|} d^4x
$$
Thus the action now depends on an unspecified function of the Ricci scalar. The case $f(R) = R$ represents classical General Relativity. It is this function that we tune in these models to get the desired properties. To test the general case rather than thest each $f(R)$ individually, you can taylor expand $f(R) = c_0 + c_1 R + c_2 R^2 + ...$ and constrain your expansion to however many degrees of freedom you think your observations can constrain. Just having $c_0$ obviously corresponds to a cosmological constant, and $c_1 = 1$ is, as we have said, the usual GR. Note however, that because this modifies gravity as a whole, and not just cosmology, other observations could be brought into the fold to constrain more parameters, and we are not limited here by the number of degrees of freedom cosmological probes can constrain. If cosmological tests of GR agree with smaller scale tests of GR in a manner consistent with an $f(R)$ modified gravity, these models could be proven to be onto something. Here of course we must look at all the tests of GR to get a fuller picture, but that is another review in and of itself. Here we'll simply note that more generally, the challenge when modifying the left hand side is to only change cosmology, but remain consistent with the smaller scale tests. The theory of general relativity is well tested on a number of scales. The accelerated expansion of the universe is the only potential anomaly when comparing general relativity to observations. To modify it such that the Friedmann equations gain an acceleration term while all other tests remain consistent is tricky business for theorists working in the field. The $\Lambda$ correction works because it essentially maintains the same theory except for adding a uniform scalar field in the theory. Hence, while its impact on the global metric is significant, it has no impact on local tests of gravity because it doesn't introduce a gradient to the energy density field. Other theories run into tuning problems precisely because they cannot naturally mimic this behavior. Chameleon theories get around this problem by adding a field that becomes extremely massive (i.e. an extremely short scale force) when surrounded by matter. Hence, its local impact on gravity is negligible in all the contexts in which general relativity has been tested. Ideally, your model would remain consistent with pre-existing tests, but still offer hope for some testable differences (ideally with a rigorous prediction set by the scale of acceleration) beyond cosmological acceleration so as to make your model more believable.
Finally, we can consider if the assumptions that went into the Friedmann equations are correct. One approach (not the one you mentioned) is to consider the impact of inhomogeneities. While it was always obvious that homogeneity did not hold below certain scales, the traditional approach has been to simply conclude that while the Friedmann equations will therefore not accurately describe the "local" scale factor, it would describe the "average" scale factor so long as we only considered it to be accurate above the inhomogeneous scales. In other words, they describe the expansion of the "average" scale factor. Mathematically speaking, this is not the rigorous way to get around this problem.
Ideally, what one would do instead, is to calculate the average of the expansion of the scale factor. Some theorists have proposed methods of doing this in ways that would add an acceleration term to the Friedmann equations. For a matter dominated universe these become :
\begin{eqnarray}
\left( \frac{\dot{a}}{a} \right)^2 &=& \frac{8 \pi G}{3} \left( \left\langle \rho \right\rangle - \frac{k}{a^2} \right) - \frac{1}{2} Q(t) \\
\frac{\ddot{a}}{a} &=& \frac{- 4 \pi G}{3} \left\langle \rho \right\rangle + Q(t) \\
\mbox{Where } Q(t) &=& \frac{2}{3} \left( \left\langle a^2 \right\rangle - \left\langle a \right\rangle^2 \right) - 2 \left\langle \sigma^2 \right\rangle
\end{eqnarray}
Where $\sigma$ depends on the level of structure formation of the universe. Because this model builds on the effect of inhomogeneities, it could explain why acceleration is a late time effect. Earlier in the universe's history, structure fromation was in its infancy and could not provide the necessary level of inhomogeneity. As the universe grows, structure formation kicks in, and so too does this effect become more significant. The biggest challenge for these models is providing a definite value for $Q$. This has proven to be a difficult task because estimating $\left\langle a^2 \right\rangle - \left\langle a \right\rangle^2$ analytically is extremely difficult, and simulating this sort of GR on such a large scale presents problems of its own.
Another potential effect of these inhomogeneities is the loss of energy as light enters large potential wells, only to exit them after they have significantly expanded. In this view, the observations would be explained by what happens to the photons along the line of sight, and no actual acceleration would take place. In some ways, this is equivalent to a sort of physically motivated tired light model. So far, these models have not been able to explain the observations assuming only a matter dominated universe. It is, however, entirely possible that a small portion of the observed effect is due to this. Keep in mind that this explanation would only affect supernovae measurements, and hence could not explain the now independent confirmation of other probes regarding acceleration. At best, they could one day reduce tension between supernovae and other probes.
The Lemaitre-Tolman-Bondi model assumes that homoegeneity is just not true at all, and that we live inside a gravitationaly collapsing system, near its center. It's a little out there, but the thing about cosmology is that it's rather easy to make far fetched models and go "hey, you never know", because hey, we don't :P It's a big universe out there and we only have access to a tiny portion of it, so why not consider I guess. The key to disproving this one will be to increase the agreement between the various probes. The metric it leads to doesn't have the same impact on all the different distance probes in the same way that $\Lambda$CDM does. So it's disprovable, and for the time being can be tuned to explain the observations, but a little far fetched to begin with. I'd guess that you can change the toy model they use to model the inhomogeneous shape of the universe and almost always get something that'll fit the data, but the more you do that the less convincing your theory becomes.
As a side note, though this isn't your main question, it's not odd that we're using an effective fluid model. It's not counting on Dark Energy having properties similar to a fluid in the traditional sense, it's just telling you that you're using certain symmetries to simplify you $T_{\mu \nu}$ tensor when solving for the Friedmann equations. Essentially, if you're in the reference frame where the momenta of the universe sum to $0$ then $T_{\mu \nu} = diag(\rho, p_x, p_y, p_z) = diag(\rho, p, p, p)$, and $w$ is just the ratio of of $\rho$ to $p$. It's just a way to simplify equations and clarify your degrees of freedom.
Best Answer
The accelerated expansion of the universe is not direct evidence for dark energy, i.e. a perfect fluid contribution to the stress-energy tensor with $w = -1$. Dark energy is just by far the simplest thing that fits the data well.
It's simple to cosmologists because they are used to dealing with matter in the form of perfect fluids, and dark energy is just another one. And it's simple to particle physicists because it can be sourced by a constant term in the Lagrangian, the simplest possible term.
At the classical level, the distinction you're making is not really important. A cosmological constant, which you call a "modification of gravity", amounts to adding a constant term to the Lagrangian. But the standard description of dark energy also amounts to adding a constant term to the Lagrangian. They're the exact same thing -- a constant is a constant, it doesn't come with a little tag saying if it's "from" gravity or something else. Neither is more elegant because functionally all the equations come out the same. It's a philosophical difference, not a real difference.
But the situation changes dramatically when you account for quantum effects. That's because we know that QFT generically produces vacuum energy, i.e. sources a constant term in the Lagrangian, whether or not we put that term in classically or not. So even if you do explain the accelerated expansion by some other mechanism, you have to explain why this one isn't in effect. This is a difficult argument to make, because the contribution from QFT is already too big even if you only trust it up to tiny energy scales like $1 \text{ eV}$!
Of course, there is room to work on alternative theories, such as quintessence and phantom energy; these are functionally different because they correspond to a perfect fluid with $w \neq -1$. At present, observational constraints show that the acceleration can only be explained by one additional perfect fluid if that component has $w = -1$ to within about 20% accuracy. If you wish, you can think of these theories as a modification of gravity by just moving their contributions to the other side of the Einstein field equation.
My impression is that these theories simply ignore the QFT vacuum energy without explanation. That's the real elephant in the room, and probably the reason so few physicists work on explaining the accelerated expansion. We have an automatic mechanism to explain the expansion, and that mechanism appears to be $10^{120}$ times more powerful than it should be. It seems premature to start to consider additional mechanisms before understanding this one first, and there seems to be no possible understanding of it now except via the much-hated anthropic principle. The ultimate explanation of the expansion will be a job for physicists of another millenium.