Frequency-dependent skin penetration depth of terahertz radiation determined by water sorption desorption

: A multilayered water–skin model is used to experimentally verify a new sensing method for determining the skin penetration depth of radiation with 0.1–0.9 terahertz (THz) frequencies. A water overlayer is dripped on a skin sample to form a multilayered structure for dynamically measuring the reflected THz-wave amplitude during water desorption. Skin penetration depths can be successfully derived by using the multilayered water–skin model and by considering the measured reflectivity, water dielectric constants, and effective thicknesses of the water overlayer on the skin sample. The maximum penetration depth is approximately 0.3 mm and is obtained with wave frequencies of 0.4–0.6 THz. This penetration depth encompasses the stratum corneum (SC) and part of the epidermis. The high penetration depth of 0.4–0.6 THz waves is also confirmed in the dried and damaged SC.


Introduction
The optical methods of microscopy and spectroscopy have been frequently adopted to noninvasively identify molecules under skin surfaces. The wavelength-dependent penetration depth of light in biological tissues has crucial applications in dermatological pharmacology and photobiological impact assessment [1]. For example, visible light is utilized in the photodynamic therapy of skin disease [2], and infrared ray (IR) is used to monitor the diffusion of drugs or the hydration level of skin [3,4].
Terahertz (THz) radiation/waves in the frequency range of 0.1-3 THz are located between microwave and far-IR regions. These waves have been widely applied in diverse fields using the mature technology of time-domain spectroscopy (TDS) [5]. THz radiation can be applied to distinguish diseased skin from normal skin based on its water-sensitive response, which appears in the time-domain waveforms or dielectric constant spectra of skin in a reflective THz-TDS [6][7][8][9]. The skin penetration depth of THz radiation determines the sensitivity of the THz-TDS system to pathological tissues. The high sensitivity of the THz-TDS system, in particular, is required for the early diagnosis of skin cancer. In studies by [10,11], the skin penetration depth of THz radiation is determined as the reciprocal of THz absorption coefficients measured from THz wave transmission in the skin. However, THz absorption coefficients that are derived from reflected THz waveforms cannot provide the exact skin penetration depth of THz radiation. The issue is that the Fresnel reflection principle of a single-layer skin surface, used to analyze the reflected waveforms, ignores wave interference in the multilayered internal structure of the skin [12]. Moreover, the frequency-dependent skin penetration depth of THz radiation has not been explored through reflective THz-TDS. Vol. 26, No. 18 | 3 Sep 2018 | OPTICS EXPRESS 22709 In this work, the frequency-dependent skin penetration depths of THz radiation are experimentally characterized by analyzing the interference effects of THz waves in a multilayered water-skin structure and are numerically verified through finite-element method (FEM) simulation. To easily observe the THz interference effect in the multilayered skin structure, a water overlayer is dripped on the surface of a porcine skin sample. When the water overlayer gradually desorbs or evaporates with time, the dynamic variations in reflected THz signals from the water-skin composite medium are measured. The frequency-dependent absorption coefficients, refractive indices, and penetration depths of THz radiation in the skin sample can be successfully derived from the measured reflectivity spectra of THz radiation on the basis of the appropriate multilayered model. In the sensing scheme, the composition of the fully hydrated internal skin structure and the water overlayer that is adsorbed on the skin surface are assumed as bulk water with known THz dielectric constants and thicknesses [13]. Comparing to the natural evaporation of the water overlayer on a normal skin, the damaged skin is obtained from the treatments of hot air and boiling water on the stratum corneum (SC) that individually modify the water content and THz dielectric properties of the skin surfaces. The reflective response of damaged skin samples to THz radiation is determined to confirm the skin penetration depths of THz radiation, which is analyzed at the discussion section of the presentation. Porcine skin is suitable for modeling human skin [14]. The skin samples used in this work were prepared from fresh skin with partial subcutaneous fat that was harvested from the back of a Yorkshire pig (female, 4 months old). Although the pig died of disease, it had healthy skin. Within 1-2 h of harvesting, the skin was cut into small samples (side length of <20 mm), sealed in polyethylene film, and stored in a -85 °C freezer. Prior to the experiment, the skin samples were thawed in a Ziplock bag for approximately 1 h. Hairs were cut as short as possible. The skin sample was then loaded in a sample holder under room temperature (28 °C) and normal atmosphere (1 atm). The weights and side lengths of the skin sample loaded in the holder are 1.0-1.4 g and 1.0-1.5 cm, respectively. Figure 1(a) illustrates the sample holder used to hold the multilayered water-skin structure. The sample holder comprises a polypropylene (PP) plastic box, a PP cap, and a sponge. When the skin sample is loaded in the PP sample holder, the sponge and cap compress the skin sample, thus preventing measurement uncertainty from irregularities and fault structures on the skin. A circular hole perforated in the center of the PP cap center is used as the operation window. In the experiment, 4.0 and 5.2 mm window diameters are applied to the normal and damaged skin samples, respectively, which are given from different conditions of surface limpness for the highest reflectivity. The skin surface obviously bulges inside the operation window to construct one well-controlled reflectance surface, thereby flattening the random wrinkle. The window diameter for the normal skin should be smaller than that of the damaged one because the tiny wrinkle of normal skin should be flattened by a sufficiently small hole, which is perforated on the PP cap. Although such the holding method is not useful for certain biotissues detected by THz waves [15], it is advantageous for THz waves to sense the skin surfaces with the optimal reflectance. Besides of the manipulation of skin surface via this operation window, water can be dripped on the sample skin surface (SC side) through it. Only the skin surface within the operation window is exposed to the ambient atmosphere and irradiated with THz waves, as illustrated in Fig. 1(b). The other section of the skin surface is isolated and covered by the cap of the sample holder to prevent evaporation.

Sample preparation
The periphery of the operation window and the joint between the cap and PP box were smeared with petrolatum. During water adsorption, the water drop that wetted and adsorbed on the skin surface formed a water overlayer within the operation window. The water overlayer can be manipulated within the operation window without leaking when measuring the THz signals reflected from the multilayered water-skin structure during water desorption. The water overlayer on or inside the skin is exhausted through natural evaporation or through exposure to hot air flow during water desorption. The temperature of the hot air flow was controlled at approximately 40 °C to completely dry the water that had adsorbed on the skin surface without causing tissue damage. Figure 1(c) presents a photograph showing that the porcine skin surface only shrinks within the operation window after hot air treatment. The circular operation window of the sample holder can precisely control the area of the water overlayer on the skin surface. The area of the water overlayer, therefore, approximates that of the operation window. The exact water overlayer area is necessary to obtain the water overlayer thickness when the weight and physical density (1.0 g/cm 3 ) of the water overlayer are used in thickness estimation.
To modify the dielectric properties of skin samples, we damaged the skin layers that include the SC and partial epidermis through treatments with boiling water, freezing at −85 °C and thawing. The SC contains a dense network of keratin, a type of protein that prevents water evaporation from the skin [16]. Keratin can be denatured through boiling, and the formation of ice crystals through freezing at −85 °C causes massive membrane disruption. Hence, the damaged skin surface after thawing drastically shrinks upon exposure to ambient atmosphere. Figure 1(d) shows that the damaged skin surface shrinks and darkens during water evaporation. By contrast, the undamaged (normal) skin exhibits completely preserved surface morphology and does not shrink during natural evaporation. Figure 2(a) schematically shows one part of the THz wave optics in a THz-TDS system, which confocally inputs and reflects one THz pulse from a skin sample or a metal mirror by using a THz plastic lens with a focal length of 25 mm. The THz pulse spectrum is shown in the Fig. 2(b) and ranged in 0.1-2.0 THz frequency. The focused beam size in the available spectral range is approximately 3 mm, and the corresponding Rayleigh length is sufficiently longer than 2.3 mm [17]. The confocal condition can maintain for the surface shrink of the skin sample because the maximum shrinkage depth of the skin surface is only 0.3 mm in prediction, which is considerably smaller than the Rayleigh length range of 2.3 mm. The metal mirror can reflect all THz waves at all powers with negligible reflection loss. Thus, the mirror-reflected power approximates the total system power to radiate the sample in the reflectivity measurement.

Sensing method
The top surfaces of the skin sample and metal mirror were irradiated with THz waves at the same z-axial focus position by adjusting the translation stage (δz) of sample holder to make their main peak positions of electric field oscillations (i.e., THz waveforms) at the same time delay. Within an equal time-domain scan range and a fixed main peak position in the THz-TDS measurement, the THz waveforms measured from the skin and metal mirror were then transformed through fast Fourier transform (FFT) to obtain the reflected power spectra of the skin and the reference mirror [18]. The FFT power spectra derived from the mirror and skin reflectors are denoted as P ref (ν) and P sample (ν), respectively. The THz reflectivity spectra of the skin, Re.(ν), can be obtained on the basis of the power ratio, i.e., P sample (ν)/P ref (ν), to illustrate the spectral properties of the target skin. During the THz wave-sensing process, one water drop of approximately 10 mg in weight was first dripped through the operation window to wet the skin surface and to form a multilayered water-skin structure via water sorption [ Fig. 1(a)]. Based on the open frame of the sample holder, the residual weight and reflected THz signal of the added water overlayer were dynamically inspected during natural evaporation. The reflected THz waveforms from the water-skin and mirror reflectors were sequentially measured at different evaporation times by exchanging the two reflectors until the water overlayer was exhausted. The thicknesses of water overlayers decrease at different evaporation times and can be used for THz wave sensing purpose.
Water overlayer thicknesses and related THz waveforms were simultaneously recorded and averaged within 12 s. Reflective THz-TDS was performed with the acquisition rate of 10 Hz for each waveform. Consequently, approximately 120 waveforms were averaged within 12 s. The total weight of the sample was also monitored in situ per 10 s to obtain the water overlayer thickness. Each reference waveform [i.e., P ref (ν)], which was obtained from the metal mirror reflection, was acquired after the sample waveform [i.e., P sample (ν)] for the estimation of THz reflectivity [Re.(ν)] from one multilayered water-skin structure. Although the Re.(ν) is not a real-time response, we assume that THz input power remains constant in the replacement action between the sample and mirror. Furthermore, the THz-TDS system enables the acquisition of the THz reflectivity spectra on the skin with a spectral range of nearly 2.0 THz [ Fig. 2(b)] in a waveform measurement without any mechanical frequency tuning. Given that the measurement criterion for a THz waveform is based on the 10 Hz scanning frequency, the water overlayer thickness used in the sensing process are constant for all THz frequencies of 0.1-2.0 THz.

Modeling
To verify the measured reflectivity on the multilayered water-skin structure during water desorption [ Fig. 1(a)], the theoretical reflectivity is calculated on the basis of the FEM results by using the COMSOL Multiphysics software. In the FEM calculation, the model of multilayered water-skin structure is illustrated in Fig. 3(a), where the thickness of an outer water overlayer is denoted as Z 2 . The skin surface that is in contact with the water overlayer is set as the origin of the z-axis. The skin thickness of the model that interacts with a THz wave is considered as one THz interactive layer and is denoted as a certain thickness of Z 3 . IR spectroscopy shows that the water content of water-rich dermis exceeds 70% [3]. Therefore, for the FEM calculation, the fully hydrated innermost skin that adjoins the interactive skin (Z 3 ) is assumed as one water-like tissue layer with an unlimited thickness Z 4 . Figure 3(b) shows a diagram for the multilayered structure of the damaged skin. There is high water occupation ratios at the outermost layer, which is much higher than that of the normal SC or epidermal layers and estimated in the discussion section of the presentation. This result corresponds with the occurrence of water sorption on skin that has been damaged through boiling water and −85 °C freezing treatments. Thus, we assume that THz wave irradiation on the surface tissue of damaged skin is equivalent to that on bulk water. As shown in Fig. 3(b), the penetration depth (Z THz depth ) at a certain THz frequency is considered as a constant during water desorption in water-rich and intrinsic layers of the damaged skin structure, i.e., Z 2 + Z 3 = constant. Hence, the intrinsic skin thickness (Z 3 ) increases, whereas the thickness of the embedded water layer (Z 2 ) decreases during evaporation. This multilayered structure of the damaged skin is used to express the interference effect of THz wave with the deepest penetration depths.   substantially higher than that of the skin under a light water overlayer (0-1.9 mg) and is nearly consistent with the theoretical Fresnel reflection of bulk water as illustrated by the black line in Fig. 4 [19]. These sensing results experimentally indicate that the reflectivity of THz waves with frequencies of 0.1-2.0 THz increases with the weight of the water overlayer. To remove the water residue adsorbed on the skin surface and to detect the THz spectral response of the shrunken skin, 4 mg of water is removed from the sample through hot air treatment. The dried skin sample is denoted as "−4 mg" in Fig. 4.

Sensing results
The response of Re. to different water amounts (W) at 0.137, 0.275, 0.413, 0.551, 0.688, and 0.826 THz is expressed by the black hollow circles in Fig. 5. The reflectivity of waves with frequencies exceeding 0.9 THz is not analyzed given the low signal-to-noise ratio in the sensing work. The maximum THz reflectivity is obtained under the outer water overlayer weight of 6-7 mg for all the observed THz frequencies in Fig. 5. When W exceeds 7 mg, Re. slightly decreases. Under a water overlayer with the weight of 0-6 mg, the Re. is approximately proportional to W, and the lowest reflectivity is obtained when the water overlayer weight is 0-2 mg. For example, 1.4 and 1.9 mg of water contribute to the lowest Re.  Using the sample preparation method (section 2.1), the associated effective thicknesses (Z eff ) of water overlayers are obtained and indicated in the upper horizontal axis versus the corresponding W values. Based on the multilayered water-skin structure [ Fig. 3(a)] and the estimated Z eff , the calculated Re. values are indicated in Fig. 5 by hollow red squares and blue triangles, which are two conditions to approach the measured results. These two conditions used in reflectivity calculation is determined from the adsorbed water amounts and explained as follows to derive the THz frequency-dependent properties in skin, i.e., absorption coefficients, refractive indices, and penetration depths. Figure 5 reveals that the distinctly high and low THz reflectivity of multilayered waterskin structure occur, respectively, for the thick and thin water overlayer. When THz waves propagate through a thick water overlayer (i.e., Z 2 thickness > 0.4 mm), the THz wave power is nearly exhausted outside the skin [ Fig. 3(a), z > 0] without obvious wave interference among Z 2 , Z 3 , and Z 4 . Thus, in the FEM calculation, the skin sample is assumed to have a certain thickness Z 3 to approach THz wave reflection with very weak interference. Such the weak interference in THz wave reflection also performs in the dry/blank skin condition, i.e., Z 2 = 0, because the reflection of skin surface is assumed not so strong to interference the secondary reflection at the inner hydrated layer Z 4 . When the water overlayer evaporates until the effective thickness of Z 2 is less than 0.2 mm but larger than 0 mm, the THz wave is transmitted through the water overlayer (Z 2 ) and penetrates the skin layer, thereby resulting in very strong wave interference among Z 2 , Z 3 , and Z 4 . To obtain the solution of the dielectric constants of skin in the multilayered water-skin structure from the reflectivity observation, one skin penetration depth should be assumed as one thickness without interference effect. Therefore, we assume a sufficiently large Z 3 value at 1 mm to derive the absorption coefficients (α eff ) and refractive indices (n eff ) of skin in theory when the outer overlayer thickness Z 2 is 0 and larger than 0.4 mm. Based on the absorption coefficients and refractive indices of skin, the skin penetration depth is then estimated using the THz reflectivity with strong wave interference in the multilayered water-skin structure [ Fig. 3(a)] when the outer overlayer thickness Z 2 is between 0 and 0.2 mm, i.e., 0 < Z 2 < 0.2 mm. In accordance with the condition of weak wave interference, the n eff and α eff of the skin sample can be estimated through the optimal approximation of the FEM-calculated reflectivity to the measured one for each water overlayer thickness Z 2 (i.e., Z eff in Fig. 5). To calculate the theoretical THz reflectivity and dielectric parameters of normal skin, the FEMbased iterative approximation method is used and initialized at the certain values of n eff and α eff [20], which are respectively shown in Figs. 6(a) and 6(b) by the black hollow circles. Black solid triangles illustrate the sensing results of n eff and α eff , respectively, in Figs. 6(a) and 6(b). The derived n eff and α eff at different THz frequencies cause the calculated reflectivity to match the measured reflectivity for the outer water overlayer thickness exceeding 0.4 mm and at 0 mm. The FEM calculated reflectivity based on the weak interference condition is illustrated by the red hollow square in Fig. 5. However, the reflectivity based on the weak interference condition of the multilayered structure obviously deviates from the measured one when the outer overlayer thickness is in the range, 0 < Z 2 < 0.2 mm, representing THz waves strongly interference for those small thicknesses of the water overlayer on the skin.
On the basis of the estimated n eff , α eff (Figs. 6(a) and 6(b)) and the layer configuration [ Fig.  3(a)] with strong interference, i.e., 0 < Z 2 < 0.2 mm, the interactive thickness (Z 3 ) of THz wave inside the skin can be evaluated via the FEM-based iterative approximation of the calculated THz reflectivity to the measured reflectivity. The interactive thickness Z 3 is obtained from two Z 2 values, 0.111 and 0.151 mm, and expressed by the black solid squares in Fig. 6(c). For the 0.688 and 0.551 THz waves, Z 3 value is not converged in the iterative approximation for the 0.151 mm of Z 2 , so that, the calculated THz reflectivity cannot be derived as shown in Fig. 5. Here, the interactive thickness Z 3 of skin model can be regarded as the skin penetration depth of THz radiation.
As shown in Fig. 5, the reflectivity under W = 0 mg is normally higher than that under W = −4 mg. This behavior is consistent with the THz response of hydrated skin [21]. We further consider the decrease in the percentages of reflectivity (δRe.) under two water desorption conditions, namely, W = −4 and 0 mg (Fig. 5) We successfully evaluate the dielectric parameters [Figs. 6(a) and 6(b)] and penetration depth Z 3 [ Fig. 6(c)] of the multilayered water-skin structures (Fig. 3) from the spectral response of THz reflectivity Re. in Fig. 5. Moreover, Re. of Fig. 5 can also reasonably explain the frequency-dependent property of skin penetration depth associated with reflection decrement δRe. from a small amount of water loss out of the skin tissue. We further compare the frequency-dependent penetration depth of skin and input THz power spectrum as expressed in Figs. 6(c) and 2(b), respectively. The largest penetration depth of 0.3 mm is obtained with 0.4-0.6 THz but not with the highest input power range of 0.6-0.8 THz. This finding implies that the spectral range of the highest THz-wave incident power on the skin does not coincide with that of the deepest penetration depth. Therefore, in a reflective THz-TDS system, the frequency-dependent property of skin penetration depth in 0.1-0.9 THz can be qualitatively interpreted as the response of THz frequency (i.e., THz photon energy) during interaction with skin tissue, instead of THz wave amplitude.

Discussion
The photonic response of skin in the sub-THz spectrum is related to the spiral ducts of sweat glands [22][23][24]. However, the photonic response of a multilayered skin structure has not been experimentally deduced from wave interference results [25]. The sensing result shown in Fig.  6(c) indicates that the frequency-dependent penetration depth of THz wave, which is correlated to wave interference, can be derived from the Re. of a multilayered water-skin structure during water desorption. In the experiment, surface-modified skin samples are additionally adopted to confirm that the maximum penetration depth is obtained under frequencies of 0.4-0.6 THz. That is, the layers of skin tissue within the penetration depth of THz wave are further modified to have different THz dielectric properties, which can be observed from the variation of reflected THz waves.
For the in situ monitoring of the THz wave response on damaged skin during water desorption, we observe the variation in the electric field amplitude instead of the power reflectivity of the waves. Figure 7(a) illustrates the average THz waveforms that are obtained at different evaporation times for durations of 10-100, 110-200, 210-300, 310-400, and 410-500 s. Within each duration, 10 waveforms are obtained and averaged as one waveform, as shown in Fig. 7(a). The sample weights are also recorded and averaged at the evaporation durations of 0-10, 100-110, 200-210, 300-310, and 400-410 s. Sample weights decrease during water evaporation. Therefore, the average and deviation values of THz wave amplitude can reasonably represent that a water layer in the skin modulates the input THz waves.   Fig. 1(d). We consider the sample weight and THz waveform measured at 0-10 and 10-100 s, respectively, as reference signals and compare them with signals measured at different evaporation times to identify the variations in different water amounts and their correlated phase delays. Figure 7(b) illustrates the optical path differences (OPDs) δ(n eff z) at different evaporation times. The OPDs are estimated from the measured phase delay shown in Fig. 7(a) on the basis of the formula δφ = 2πνδ(n eff z)/C, where z, n eff. , ν, and C are the wave propagation distance along the z-axis [ Fig. 2(a)], the effective refractive index of the skin, the frequency of the THz wave, and the speed of light in the air, respectively. The physical propagation path of the THz wave (z) is extended from the water level in the skin. The shrunken skin depth, estimated from the OPD value, is larger than the effective thickness of the embedded water layer, derived from the water weight loss within the operation window of 5.2 mm diameter (section 2.1). The effective thicknesses of the water layers embedded in the damaged skin surface [ Fig. 3(b), Z 2 ] are 0.012, 0.042, 0.047, and 0.078 mm for the water weights of 0.33, 1.12, 1.25, and 2.08 mg [ΔW, Fig. 7(b)], respectively, and individually present shrunken depth [i.e., δ (n eff z), OPD] of 55%, 93%, 75%, and 80%. Other fractional thickness of the shrunken skin is considered as the intrinsic skin tissue [ Fig. 3(b)]. Such the high water percentages within the skin outermost layer (i.e., shrink depth or OPD) experimentally validate one water-rich layer at the outermost skin [ Fig. 3(b)] to explain the deepest penetration depths in Fig. 6(c). Figure 8 illustrates the variation in the electric field amplitude (ΔE) of damaged skin with different water losses. ΔE values are obtained from the THz-FFT spectrum shown in Fig. 7(a) and are compared with the reference spectra of the averaged waveform obtained during 10-100 s of evaporation [ Fig. 7(a), black line]. The response between ΔE and ΔW (or Z 2 ) at different THz frequencies is based on the linear fit results (Fig. 8), where the parameter of the adjusted determination coefficient R 2 (i.e., the adjusted R 2 or adj. R 2 ) shows how well the data points fit the line. The adjusted R 2 of the linear fit thus evaluates the proportional relation between ΔE and ΔW. Comparing to R 2 value, the adjusted R 2 value is advantageous to analyze small amounts of data points without bias on the determination coefficient of the linear fit. In the case (Fig. 8), there are only four points at each frequency and thus the adjusted R 2 values are used in the statistical analysis of THz wave response in the desorption process of a damaged skin sample. A high value of adjusted R 2 indicates that the ΔE value is stably changed by a water amount ΔW under water evaporation. By contrast, a low adjusted R 2 indicates that ΔE easily fluctuates for the same water amount. Figure 8 illustrates that the waves with frequencies of 0.383, 0.460, and 0.536 THz have a considerably high adjusted R 2 value of greater than 0.9. The THz waves with the high adjusted R 2 also have distinctly high penetration depths of 0.27-0.31 mm, as shown in Fig. 6(c). For other THz waves with penetration depths of less than 0.27 mm (>0.6 THz or <0.38 THz), the adjusted R 2 decreases under the same ΔW condition. Therefore, the experiment to detect the water layer under the damaged skin surface using 0.1-0.9 THz waves suggests that penetration depths are correlated to THz frequency and that the maximum penetration depth is obtained with 0.4-0.6 THz frequency.
Given that ΔE values are dynamically measured and averaged over 90 s of evaporation, the error bars (i.e., ΔE variation range) in Fig. 8 are apparent and different at various THz frequency due to the emission power of THz system (Fig. 2). For example, THz waves nearby the spectral peak, 0.7 THz [ Fig. 2(b)], have the smallest ΔE error bars, where the 0.28-0.53 unit of error bar is found in the frequency range of 0.690-0.843 THz (Fig. 8). THz waves in the spectral edge, 0.230-0.383 THz [ Fig. 2(b)], have the considerably large ΔE error bars, 11.2-25.4 unit (Fig. 8). It represents that high-power THz waves perform stable reflection and the low power waves are oppositely not stable with large fluctuation. The error bar size, detected by one THz pulse, is not random among the THz frequencies because no additional deviation source occurs during the frequency adjustment. It thereby provides the background of frequency-dependent power stability during the 90 s observation and can be reduced as small as possible in a short term observation. Each ΔE value is averaged by 10 waveforms within the long term observation (90 s) and each waveform is also averaged by the timedomain scanning frequency about 10 Hz. One ΔE value in Fig. 8 corresponds to be averaged by 100 waveforms during the 90 s, and reliable to respond ΔW on the damaged skin for the frequency-dependent R 2 values and slopes, which is not dominated by the power stability background, i.e., ΔE variation. The ΔE response at 0.383 and 0.460 THz is proportional to the amount of water loss. This relationship is consistent with the δRe. response of normal skin after the removal of approximately 4 mg of water as shown in Figs. 5 and 6(c). That is, the water content of skin proportionally increases the reflectivity of THz waves. By contrast, the ΔE response at 0.536 THz shows that high THz reflectivity is associated with low water content (middle figure in Fig. 8). The inversely proportional response of ΔE vs. ΔW results from the constructive interference of THz radiation in damaged tissue. Figure 9(a) schematically illustrates the interference effect between the first and second reflections of the input THz wave when the embedded water is nearly exhausted through evaporation or when Z 2 [Figs. 3(b) and 9(a)] is approximately 0 mm. While the first and second reflected waves are in phase at the air space, the criterion of the constructive interference is the OPD inside skin approximates one wavelength. When the water loss (ΔW) increases through evaporation, the thickness of the intrinsic skin (Z 3 ) increases [Figs. 3(b) and 9(a)] to match the wavelength of 0.536 THz, thus exerting constructive interference at the system detector. Figure 9(b) shows the FEM-calculated reflectivity values (Re.) of waves with frequencies of 0.460 and 0.536 THz under different water thicknesses (Z 2 ) within a fixed penetration depth [Z THz depth , Fig. 3(b)]. Several Z THz depth values of the damaged skin are selected for numerical simulation. These values include 0.15, 0.25, 0.32, 0.40, and 0.45 mm [ Fig. 9(b)]. In the FEM calculation, the refractive indices and absorption coefficients of the bulk water and skin layers are, respectively, 2.35, 1.5 RIU and 165, 21 cm −1 [19]. Given that the fitting lines of 0.460 and 0.536 THz frequencies in Fig. 8 intersect at -0.02-0 mm-Z 2 , we consider it as a recognized feature in Fig. 9(b) to obtain the 0.32-0.40 mm-Z THz depth as the solution range in FEM calculation. This FEM calculation result validates that the maximum skin penetration depth (Z THz depth ) of more than 0.3 mm, but not exceeding 0.45 mm, is achieved by waves with frequencies of 0.4-0.6 THz.

Conclusion
The THz frequency-dependent property of skin penetration depth is experimentally identified during water sorption-desorption on a skin sample. For water sorption, the multilayered water-skin medium is constructed to observe THz wave interference effect by placing a water drop on the sample skin surface. Water desorption is considered as the natural evaporation of water under ambient atmosphere until the skin surface is dry. The multilayered water-skin structure is continuously irradiated with THz waves during water evaporation to dynamically record the reflected amplitudes of THz waves. Skin samples under water overlayer thicknesses of more than 0.4 mm and the dry condition can be modeled as a multiple layer configuration with weak interference of THz wave. The THz refractive indices and absorption coefficients of the weak interference condition can thus be derived via the FEM-based iterative approximation of THz reflectivity. Skin under water overlayer thicknesses of less than 0.2 mm but larger than 0 mm can be modeled as a multiple layer configuration with strong interference of THz wave. The THz interactive layer thicknesses correspond to THz penetration depth and are consequently obtained through the iterative approximation of FEM calculation of THz reflectivity. The sensing results show that the 0.1-0.9 THz waves have skin penetration depths of 0.1-0.3 mm and those waves with 0.4-0.6 THz frequencies especially have the maximum skin penetration depth of 0.3 mm. To confirm the maximum penetration depths of THz waves and to validate the multilayered skin model, the skin surface is further damaged with boiling water and freezing at -85 °C to induce massive membrane disruption, thereby forming a large porous space in the skin surface with a high water content. The porosity or water content of the damaged skin tissue with a thickness of 0.1 mm is nearly 80% and thus approximates a water-like tissue layer at the skin surface. When THz waves reflect from the damaged skin, the waves with frequencies of >0.6 THz or <0.38 THz possibly have penetration depths of less than 0.27 mm and could be blocked by the outermost damaged tissue, thereby performing high reflected field fluctuation during the water loss. Contrarily, the reflected field variation of THz waves with frequencies at 0.383, 0.460, and 0.536 THz are linearly related with water loss, which probably results from the penetration depths of more than 0.27 mm. In contrast to that of 0.383 and 0.460 THz waves, the amplitude variation of 0.536 THz wave is inversely proportional to the water loss amounts given the constructive interference between the first and second reflections at the outermost water layer and innermost hydrated skin tissue, respectively.