  In this letter, we provide counter-examples to demonstrate that it is
possible to retain second-order accuracy using Sakurai et al.'s method, even
when different flux boundary conditions are imposed on multiple interfaces that
do not conform to the Cartesian grid. We consider both continuous and
discontinuous indicator functions in our test problems. Both indicator
functions yield a similar convergence rate for the problems considered here. We
also find that the order of accuracy results for some of the cases presented in
Sakurai et al. are not reproducible. This is demonstrated by re-considering the
same one- and two-dimensional Poisson problems solved in Sakurai et al. in this
letter. The results shown in this letter demonstrate that the spatial order of
accuracy of the flux-based VP approach of Sakurai et al. is between
$\mathcal{O}$(1) and $\mathcal{O}$(2), and it depends on the underlying
problem/model. The spatial order of accuracy cannot simply be deduced a priori
based on the imposed flux values, shapes, or grid-conformity of the interfaces,
as concluded in Sakurai et al. Further analysis is required to understand the
spatial convergence rate of the flux-based VP method.
