> so far, I can't see any cells being involved when it comes to this method

By "this method" do you mean `surfaceDensities()`?  True, there is no `cells` input parameter (or any other, for that matter) to that method.  That's "simply" because we always assumed that there would be a single PVT region so that there would be no need for spatial variability in that method.

Your proposal, on the other hand, is as far as I can tell quite a bit more wide-reaching.  Please correct me if I've misunderstood, but aren't you proposing that all property methods accept region indices rather than cells?  That, in particular, means that the `viscosity()` method change from
```C++
void viscosity(int n, const double* p, const double* z, const int* cells,
               double* mu, double* dmudp);
```
to
```C++
void viscosity(const int region, const double p, const double* z,
               double& mu, double* dmudp);
```
and that call sites be updated accordingly (in particular, that there be separate calls for each cell).  Moreover, that's nowhere near sufficiently complete because it does not allow for end-point scaling or hysteresis effects that can (and frequently do) vary *a lot* more than per region.

You should also note that even if you *do* pass the region index, you must *still* assume that the property is per-cell.  There is, for instance, no generally good choice for the region index in the case of an interface that happens to separate two distinct regions.