## Materials and Methods
### Governing Equations  
The complete two-layer skin model contained the epidermis, the dermis, and the blood vessels (Figure 1). A 1.5 mm \(\times\) 1.5 mm \(\times\) 1.0 mm computational domain of the skin with a reconstructed blood vessel buried at a depth of 0.5 mm was used. The thicknesses of the epidermis and the dermis were set respectively as \(H_{e}\) = 60 \(\upmu\)m and \(H_{d}\) = 940 \(\upmu\)m based on statistical average values [15]. As mentioned in Section 2.1, in validating the damage model, the reconstructed blood vessels in accordance with the DSC model were embedded in the skin without the epidermis.  
In the present work, light propagation and energy deposition in the skin tissue were simulated by using our in-house TMC code [28] through which the body-fitted material grid of the tetrahedra was generated by assigning a value representing a tissue type. No interface tetrahedral cells existed, thereby avoiding the error of photon reflection/refraction in the voxel-based Monte Carlo method [29]. Photon propagations in TMC were expressed as probability distributions that describe the step size of photon movement between sites of photon-tissue interaction. Once a photon has reached an interaction site, a fraction of the photon weight was absorbed in the local tetrahedron. The \"random walk\" of photon could be simulated repeatedly until a sufficient quantity of photons was absorbed in tetrahedral cells to establish reliable trends of energy deposition. In this work, the laser beam was a flat top profile with a diameter of 1.0 mm. \(10^{8}\) photons were launched at \(z=0\)\(\upmu\)m with their \(x\) and \(y\) determined by the random number restricted in the range of the beam. The directional cosines (\(\mu_{x}\), \(\mu_{y}\), \(\mu_{z}\)) were initialized to (0, 0, 1) and the photon weight is initialized to 1.  
Energy deposition in different cells determines the heat source and thus strongly influences mass and heat flow in tissues. The following numerical formulation was implemented in the finite-volume-based CFD package FLUENT 17.0 (Ansys Inc., Canonsburg, U.S.A.) to simulate blood flow and heat transfer in vessels. The tetrahedral grid used in CFD computation was the same as that in TMC, and the energy deposition computed by TMC was treated as heat source to the governing equations. Considering that blood is an incompressible Newtonian fluid, we can write a single set of governing equations, including continuity, momentum, and energy equations in the portion of the vessel, as follows:  
\[\frac{\partial\rho}{\partial t}+\
abla\cdot(\rho\mathbf{u})=0, \tag{1}\]  
\[\frac{\partial(\rho\mathbf{u})}{\partial t}+\
abla\cdot(\rho\mathbf{u}\mathbf{ u})=-\
abla p+\
abla\cdot\left[\mu\left(\
abla\mathbf{u}+\
abla\mathbf{u}^{T} \right)\right]+\rho\mathbf{g}+\mathbf{F}, \tag{2}\]  
\[\frac{\partial(\rho c_{p}T)}{\partial t}+\
abla\cdot(\rho c_{p}\mathbf{u}\;T) =\
abla\cdot(k\;\
abla T)+S_{L}+S_{E}, \tag{3}\]  
Figure 1: Computational domain of the two-layer skin model with a reconstructed blood vessel.  
where \(\mathbf{F}\) is the surface tension, \(S_{L}\) is the energy source computed by TMC through laser pulse, and \(S_{E}\) is the energy source term related to evaporation or condensation.  
By contrast, the epidermis and the dermis were treated as a solid phase, and only heat transfer caused by laser heating was involved:  
\[\frac{\partial(\rho c_{p}T)}{\partial t}+\
abla\cdot(\rho c_{p}\mathbf{u}\ T)= \
abla\cdot(k\
abla T)+S_{L}. \tag{4}\]  
The liquid-vapor mass transfer caused by laser heating (evaporation and condensation) is governed by the vapor transport equation [7], and the inbuilt VOF module in FLUENT 17.0 was adopted to capture the interface when bubbles formed in the blood.