One more perspective: after performing a diffeomorphism that locally flattens the boundary, we get an equation of the form $$a_{ij}(x)w_{ij} + b_i(x)w_i = f(x_n) \text{ in } B_1,$$ where the coefficients are smooth (with bounds on derivatives depending only on $\Omega$) and $f$ is the Heaviside step function. Assume first for simplicity that $f$ is smooth (so that $w$ is smooth) and bounded between $0$ and $1$. Calderon-Zygmund theory gives $w \in W^{2,\,p}$ for any $p$ with an estimate. Differentiating the equation in the $e_k$ direction for $k < n$ we see that $a_{ij}(w_k)_{ij}$ is a linear combination of $D^2w, \nabla w$. Calderon-Zygmund theory thus gives $W^{2,p}$ estimates for the horizontal derivatives of $w$ for any $p$. One can continue to differentiate horizontally and apply the Calderon-Zygmund estimates to get bounds on the derivatives of $w$ of all orders in the horizontal directions, independent of the regularity of $f$. (One can equally well use basic energy estimates in $L^2$ to get the same result, if one wants to avoid Calderon-Zygmund theory). By approximation, the same holds when $f$ is the heaviside function. Reversing the diffeomorphism we see that $V$ is smooth along the boundary.

One more perspective: after performing a diffeomorphism that locally flattens the boundary, we get an equation of the form $$a_{ij}(x)w_{ij} + b_i(x)w_i = f(x_n) \text{ in } B_1,$$ where the coefficients are smooth (with bounds on derivatives depending only on $\Omega$) and $f$ is the Heaviside step function. Assume first for simplicity that $f$ is smooth (so that $w$ is smooth) and bounded between $0$ and $1$. Calderon-Zygmund theory gives $w \in W^{2,\,p}$ for any $p$ with an estimate. Differentiating the equation in the $e_k$ direction for $k < n$ we see that $a_{ij}(w_k)_{ij}$ is a linear combination of $D^2w, \nabla w$. Calderon-Zygmund theory thus gives $W^{2,p}$ estimates for the horizontal derivatives of $w$ for any $p$. One can continue to differentiate horizontally and apply the Calderon-Zygmund estimates to get bounds on the derivatives of $w$ of all orders in the horizontal directions, independent of the regularity of $f$. By approximation, the same holds when $f$ is the heaviside function. Reversing the diffeomorphism we see that $V$ is smooth along the boundary.

One more perspective: after performing a diffeomorphism that locally flattens the boundary, we get an equation of the form $$a_{ij}(x)w_{ij} + b_i(x)w_i = f(x_n) \text{ in } B_1,$$ where the coefficients are smooth (with bounds on derivatives depending only on $\Omega$) and $f$ is the Heaviside step function. Assume first for simplicity that $f$ is smooth (so that $w$ is smooth) and bounded between $0$ and $1$. Calderon-Zygmund theory gives $w \in W^{2,\,p}$ for any $p$ with an estimate. Differentiating the equation in the $e_k$ direction for $k < n$ we see that $a_{ij}(w_k)_{ij}$ is a linear combination of $D^2w, \nabla w$. Calderon-Zygmund theory thus gives $W^{2,p}$ estimates for the horizontal derivatives of $w$ for any $p$. One can continue to differentiate horizontally and apply the Calderon-Zygmund estimates to get bounds on the derivatives of $w$ of all orders in the horizontal directions, independent of the regularity of $f$. (One can equally well use basic energy estimates in $L^2$ to get the same result, if one wants to avoid Calderon-Zygmund theory). Taking a limit, we see that the same holds when $f$ is the heaviside function. Reversing the diffeomorphism we see that $V$ is smooth along the boundary.

One more perspective: after performing a diffeomorphism that locally flattens the boundary, we get an equation of the form $$a_{ij}(x)w_{ij} + b_i(x)w_i = f(x_n) \text{ in } B_1,$$ where the coefficients are smooth (with bounds on derivatives depending only on $\Omega$) and $f$ is the Heaviside step function. Assume first for simplicity that $f$ is smooth (so that $w$ is smooth) and bounded between $0$ and $1$. Calderon-Zygmund theory gives $w \in W^{2,\,p}$ for any $p$ with an estimate. Differentiating the equation in the $e_k$ direction for $k < n$ we see that $a_{ij}(w_k)_{ij}$ is a linear combination of $D^2w, \nabla w$. Calderon-Zygmund theory thus gives $W^{2,p}$ estimates for the horizontal derivatives of $w$ for any $p$. One can continue to differentiate horizontally and apply the Calderon-Zygmund estimates to get bounds on the derivatives of $w$ of all orders in the horizontal directions, independent of the regularity of $f$. (One can equally well use basic energy estimates in $L^2$ to get the same result, if one wants to avoid Calderon-Zygmund theory). By approximation, the same holds when $f$ is the heaviside function. Reversing the diffeomorphism we see that $V$ is smooth along the boundary.