A Method to Determine Human Skin Heat Capacity Using a Non-Invasive Calorimetric Sensor

A calorimetric sensor has been designed to measure the heat flow dissipated by a 2 × 2 cm2 skin surface. In this work, a non-invasive method is proposed to determine the heat capacity and thermal conductance of the area of skin where the measurement is made. The method consists of programming a linear variation of the temperature of the sensor thermostat during its application to the skin. The sensor is modelled as a two-inputs and two-outputs system. The inputs are (1) the power dissipated by the skin and transmitted by conduction to the sensor, and (2) the power dissipated in the sensor thermostat to maintain the programmed temperature. The outputs are (1) the calorimetric signal and (2) the thermostat temperature. The proposed method consists of a sensor modelling that allows the heat capacity of the element where dissipation takes place (the skin) to be identified, and the transfer functions (TF) that link the inputs and outputs are constructed from its value. These TFs allow the determination of the heat flow dissipated by the surface of the human body as a function of the temperature of the sensor thermostat. Furthermore, as this variation in heat flow is linear, we define and determine an equivalent thermal resistance of the skin in the measured area. The method is validated with a simulation and with experimental measurements on the surface of the human body.


Introduction
The study of the thermal dissipation of the human body is of great interest in multiple fields. In air conditioning projects, it is necessary to know the dissipation of the occupants according to their activity. In the study of a subject's metabolism, thermal dissipation is determined indirectly, measuring the absorbed VO2 or VCO2 emitted by the subject [1]. In the field of human physiology, all available data and techniques are used. Contact and remote thermometry are irreplaceable tools. There are numerous publications that provide temperature data of different areas and organs of the human body [2]. The thermal properties (thermal conductivity and specific heat capacity) of different parts of the human body are also of interest. Heat capacity is mainly determined from tissue sample analysis using techniques such as differential scanning calorimetry (DSC). These results are used for the study of pathologies [3,4] and the development of new bio-thermal models [5][6][7].

The Calorimetric Sensor
The experimental system has already been described in previous works [19,20], but in order to explain the calorimetric model proposed in this work, it is necessary to briefly describe the elements that constitute the sensor. The core of the sensor is a measurement thermopile (part 2 in Figure 1) located between an aluminium plate (part 1 in Figure 1) and a thermostat (part 3 in Figure 1). The thermostat is a 10 x 10 x 4 mm 3 aluminium block that has an RTD (resistance temperature detector) and a heating resistance inside. The power dissipated in this resistance is determined by a PID controller (proportional-integral-derivative controller) that maintains the programmed temperature.
The thermostat has a cooling system based on a Peltier effect cooling thermopile, a heatsink and a fan (parts 4, 5 and 6 in Figure 1). The sensor assembly, excluding the heatsink-fan block, has a thermal insulation element (part 7 in Figure 1) to decrease external disturbances.
The purpose of the sensor is to measure the heat flux (W1) transmitted by conduction from the surface of the human body to the thermostat. As a consequence of this heat flow and due to the Sebeeck effect, the measurement thermopile provides the calorimetric signal (y). On the other hand, this heat flow will produce a modification of the power dissipated in the heating resistance (W2) that keeps the temperature of the thermostat under control (T2).
To calibrate the sensor, a base has been constructed consisting of a block of insulating material (part 8 in Figure 1) that has a copper plate that contains a resistance for Joule calibration (part 9 in Figure 1). In addition, the base has a magnetic clamping system (part 10 in Figure 1) that facilitates sensor handling. Sensor calibration is performed by dissipating known powers at the base while the temperature of the sensor thermostat is controlled. Calibration measurements are designed trying to reproduce the measurement conditions in the human body. The recalibration of the sensor is completed when the sensor is disassembled for technical reasons or when it is necessary to verify its operation.  (6) fan, (7) thermal insulation, (8) calibration base (insulating material), (9) copper plate containing the calibration resistance, (10) magnets to hold the sensor in the base. The thermostat has a cooling system based on a Peltier effect cooling thermopile, a heatsink and a fan (parts 4, 5 and 6 in Figure 1). The sensor assembly, excluding the heatsink-fan block, has a thermal insulation element (part 7 in Figure 1) to decrease external disturbances.
The purpose of the sensor is to measure the heat flux (W 1 ) transmitted by conduction from the surface of the human body to the thermostat. As a consequence of this heat flow and due to the Sebeeck effect, the measurement thermopile provides the calorimetric signal (y). On the other hand, this heat flow will produce a modification of the power dissipated in the heating resistance (W 2 ) that keeps the temperature of the thermostat under control (T 2 ). To calibrate the sensor, a base has been constructed consisting of a block of insulating material (part 8 in Figure 1) that has a copper plate that contains a resistance for Joule calibration (part 9 in Figure 1). In addition, the base has a magnetic clamping system (part 10 in Figure 1) that facilitates sensor handling. Sensor calibration is performed by dissipating known powers at the base while the temperature of the sensor thermostat is controlled. Calibration measurements are designed trying to reproduce the measurement conditions in the human body. The recalibration of the sensor is completed when the sensor is disassembled for technical reasons or when it is necessary to verify its operation.
The method of measurement in the human body is as follows: initially the sensor is located in the calibration base until the thermal equilibrium is reached for the programmed temperature of the thermostat (initial steady state). The time required is 10 min from when the system is turned on. It is then placed on the surface of the skin where the measure will take place. The time applied to the skin depends on the type of measurement to be performed. A measurement with the constant thermostat temperature to determine only the heat flow requires 5 min, but a measurement with a variable temperature of the thermostat to determine, in addition, the thermal properties of the skin, requires 12.5 min. Finally, the sensor is placed back in the calibration base.
There are clear differences between calibration measurements and measurements on human skin, since in the first case the sensor is permanently located at the base, and in the measurements on human skin there is some movement of the sensor from the base to the surface of the human body. Such differences are important and will be tackled throughout this work.

Calorimetric Model
The device is modelled as a two-inputs, two-outputs system. The inputs are the power that passes through the sensor and that is to be measured (W 1 ), and the power dissipated in the thermostat (W 2 ) that maintains the programmed temperature of the thermostat. The outputs are the calorimetric signal (y) and the thermostat temperature (T 2 ).
The signal-to-noise ratio of the output signals in an experimental measurement on human skin varies from 300 to 3000, (50 dB to 70 dB). This signal-to-noise ratio allows the identification of up to two poles of the transfer functions (TFs) that link the inputs to the outputs. For this reason, a simple two-element model is proposed (Figure 2). This modelling, called localized constants, is frequently used in heat conduction calorimeters [30][31][32]. The first element represents the domain where the heat dissipation to be measured takes place (W 1 ). When the sensor is applied to the human body, this domain is the human skin. When the sensor is in the calibration base, this domain is the copper plate where the Joule calibration resistor is located. The second element represents the thermostat. The power dissipated in this domain (W 2 ) is determined by a PID controller that holds the programmed temperature. The elements have C 1 and C 2 heat capacities. As the thermal conductivity of these elements is infinite, this implies considering a temperature in each domain (T 1 and T 2 ). These elements are thermally coupled to the outside and to each other by P 12 , P 1 and P 2 thermal conductance couplings. The outside is configured by the environment that is at a T room temperature, and the cooling system that is at a T cold temperature (see Figure 2).
The power developed in each element (W 1 and W 2 ) is equal to the accumulated heat power in the domain plus the power transmitted by conduction to the neighbouring domains. The following equations (Equation (1)) describe the energy balance that relates the powers to the temperatures of each domain [19].
The calorimetric signal is a linear combination of the temperatures of the bodies that are in contact with each face of the measuring thermopile. That is, y = k(T 1 − T 2 ).
Sensors 2020, 20, 3431 5 of 20 power dissipated in this domain (W2) is determined by a PID controller that holds the programmed temperature. The elements have C1 and C2 heat capacities. As the thermal conductivity of these elements is infinite, this implies considering a temperature in each domain (T1 and T2). These elements are thermally coupled to the outside and to each other by P12, P1 and P2 thermal conductance couplings. The outside is configured by the environment that is at a Troom temperature, and the cooling system that is at a Tcold temperature (see Figure 2).  Considering that the ambient temperature and the temperature of the cold focus (T room and T cold ) are constant during the measurement and correcting the baselines of the variables (T 1 , T 2 , W 1 and W 2 ), we obtain an equivalent system of equations (Equation (2)) in which T room and T cold do not appear. In this case the values of the curves ∆T 1 , ∆T 2 , ∆y, ∆W 1 and ∆W 2 are zero at the initial instant [19].
When the sensor is placed on its base, the ambient temperature T room remains practically constant. However, when the sensor is applied to the skin there are two significant changes. The first one refers to C 1 ; now, this element is not the base but the area of the skin where the dissipation is being measured. This heat capacity is unknown and different in each area of the human body. The second change is related to the ambient room temperature around the sensor. This temperature in the neighbourhood of the skin is higher than the ambient temperature in the environment of the base. When the baselines of the experimental curves are corrected, the initial situation is corrected, eliminating the terms T room and T cold at the base, but the "new ambient temperature" in the neighbourhood of human skin is not taken into account. Therefore, we assume that the ambient temperature in the vicinity of the skin is higher by a ∆T 0 value.
We assume that this increase in temperature ∆T 0 also affects the temperature of the cold focus, since T cold = T room + ∆T Peltier , where ∆T Peltier is the decrease in temperature produced by the cooling system. The value of ∆T Peltier depends on the supply voltage (V Peltier ) of the thermopile and comes from Equation (3), whose parameters are determined experimentally [19,21].
In our case, a 0 =1.7 K, a 1 = −9.2 KV −1 ; a 2 = 1.0 KV −2 . In each measure V Peltier is constant and consequently ∆T Peltier is also constant. Thus, if there is a change in T cold , it is due only to changes in T room . When the sensor is placed on the skin, the outer surface of the thermal insulation maintains the initial temperature and gradually adapts to the new situation. This thermal adaptation of the surface is not instantaneous, and in the calculation, we will assume that this adaptation follows a variation given by  (4)). The validity of this expression is checked by adjusting the experimental curves with the curves determined by the model.
The parameters A, B and the time constant τ T are determined in each measurement. Time (t) begins when the sensor is applied to the skin and ends (t end ) when the sensor is placed back on its base. The new equations (Equation (5)) incorporate this variation in the ambient temperature (∆T 0 ).
By sorting this system of equations, and applying and operating the Laplace transform, six TFs (Equation (6)) are obtained that relate the inputs ∆W 1 (s), ∆W 2 (s) and ∆T 0 (s) with the outputs ∆Y(s) and ∆T 2 (s), where s is the Laplace variable. It is assumed that the initial value of all variables is zero and its derivative at the time origin is also zero and, for this reason, it is important to have a good initial baseline.
As can be seen in (Equation (6)), the denominators of the six TFs are the same, which implies that the poles of all the TFs are the same. However, the numerators of each TF are different. It can be verified that the transient response of the calorimetric signal and the temperature of the thermostat depend on the heat capacity C 1 where dissipation occurs. In general, an increase in C 1 implies higher values of the time constants, which results in a slower transient response. Nevertheless, this will be discussed further in the simulation section.

Identification of the Model
The hypothesis proposed in the model implies that the parameters C 2 , P 1 , P 12 , P 2 and k are characteristic of the sensor and invariant. On the other hand, the calorimetric signal (y), the temperature of the thermostat (T 2 ) and the power dissipated in the thermostat (W 2 ) are known. However, C 1 depends on the area where dissipation W 1 takes place. This power W 1 is only known when performing calibration measurements: these measurements allow determining the invariant parameters of the model.
For calibration, the sensor is placed on its base and measurements similar to those made on human skin are programmed. Initially, a constant temperature of the thermostat T 20 is programmed, and when the steady state is reached, a constant power W 10 in the base is dissipated for 300 s. Then, for 150 s, a linear increase in the thermostat temperature is programmed until T 2end and simultaneously the dissipation in the base is linearly decreased from W 10 to W 1end . Finally, and for 300 s, the temperature of the thermostat T 2end is kept constant, also keeping the dissipation W 1end constant.
In Figure 3 the experimental curves of the input powers (W 1 and W 2 ) and the experimental curves of the output signals are represented: the calorimetric signal (y) and the thermostat temperature (T 2 ). Two cases are shown: in the first one (green curves) the temperature increase in the thermostat is 10 K (from 26 to 36 • C); and in the second case (blue curves) the increase is 5 K (from 26 to 31 • C). In both cases the power W 1 dissipated at the base (red curve) is 300 mW at the beginning and 100 mW at the end.  A minimization method based on the Nelder-Mead simplex algorithm [33][34][35] is used to identify the model parameters. The error criterion used is the sum of the roots of the root mean square errors (RMSE), weighted (Equation (7)). The root mean square error is calculated with the experimental curves (subscript exp) and those calculated by the model (subscript cal).
N is the number of points and  is a weighting coefficient (N = 1800,  = 4). The calculated curves are determined with the system of equations (Equation (6)). The sampling period is 0.5 s.
The calculation starts with some initial values of the model parameters (C1, C2, P1, P2, P12 and k). With these parameters the TFs that relate the inputs to the outputs (Equation (6)) are constructed, and the outputs calculated with the model (ycal and T2cal) are determined for the known inputs W1 and W2. Then the error criterion is applied (Equation (7)) and the Nelder-Mead method that provides new parameter values (C1, C2, P1, P2, P12 and k) is applied. The iterative process decreases the error (Equation (7)) until the fit between the calculated and experimental curves cannot be further improved. Figure 4 shows a diagram of the calculation procedure. A minimization method based on the Nelder-Mead simplex algorithm [33][34][35] is used to identify the model parameters. The error criterion used is the sum of the roots of the root mean square errors (RMSE), weighted (Equation (7)). The root mean square error is calculated with the experimental curves (subscript exp) and those calculated by the model (subscript cal).
N is the number of points and α is a weighting coefficient (N = 1800, α = 4). The calculated curves are determined with the system of equations (Equation (6)). The sampling period is 0.5 s.
The calculation starts with some initial values of the model parameters (C 1 , C 2 , P 1 , P 2 , P 12 and k). With these parameters the TFs that relate the inputs to the outputs (Equation (6)) are constructed, and the outputs calculated with the model (y cal and T 2cal ) are determined for the known inputs W 1 and W 2 . Then the error criterion is applied (Equation (7)) and the Nelder-Mead method that provides new parameter values (C 1 , C 2 , P 1 , P 2 , P 12 and k) is applied. The iterative process decreases the error (Equation (7)) until the fit between the calculated and experimental curves cannot be further improved. Figure 4 shows a diagram of the calculation procedure. Table 1 shows the values of the parameters (C 1 , C 2 , P 1 , P 2 , P 12 , k) obtained in the calibration measurements. These results correspond to a series of 20 measurements in which the power dissipated in the base (400, 300 and 100 mW) and the increase in temperature of the thermostat (5 and 10 K) have been varied. The maximum RMSE in the calorimetric signal fittings is 14 µV over a 100mV peak-to-peak signal; and 4 mK on a 10 K peak-to-peak signal at the thermostat temperature fitting.   Table 1 shows the values of the parameters (C1, C2, P1, P2, P12, k) obtained in the calibration measurements. These results correspond to a series of 20 measurements in which the power dissipated in the base (400, 300 and 100 mW) and the increase in temperature of the thermostat (5 and 10 K) have been varied. The maximum RMSE in the calorimetric signal fittings is 14 µV over a 100mV peak-to-peak signal; and 4 mK on a 10 K peak-to-peak signal at the thermostat temperature fitting.  (7))

Simulation.
Calibration of the sensor allows determining the model parameters that best fit the experimental measurements (Table 1). In this section the operation of the sensor is studied from simulations. (Equation (8)) shows the relationship between the inputs and outputs of the model by means of six transfer functions (TFi) whose expressions are indicated in the system of equations in (Equation (6)). In the calorimetric response, we observe that when C1 increases, the system responds slower, which is verified in TF1 and TF3; however, the dynamics of TF2 remains unchanged. On the other hand, the dynamic response of the thermostat temperature is less affected with the variation of C1, but as expected, the response of TF4 is slower with increasing values of C1. TF5 and TF6 are less altered with the variation of C1.  (5) and (6)) obtained in the calibration measurements. Parameter

Simulation
Calibration of the sensor allows determining the model parameters that best fit the experimental measurements (Table 1). In this section the operation of the sensor is studied from simulations. (Equation (8)) shows the relationship between the inputs and outputs of the model by means of six transfer functions (TF i ) whose expressions are indicated in the system of equations in (Equation (6)). In the calorimetric response, we observe that when C 1 increases, the system responds slower, which is verified in TF 1 and TF 3 ; however, the dynamics of TF 2 remains unchanged. On the other hand, the dynamic response of the thermostat temperature is less affected with the variation of C 1 , but as expected, the response of TF 4 is slower with increasing values of C 1 . TF 5 and TF 6 are less altered with the variation of C 1 .

Simulations in the Calibration Base and in the Human Body
The simulations are performed from the systems of equations (Equation (5) or Equation (6)). First, we simulate the operation of the sensor when it is located in its base, which is completed for different heat capacities of the base: C1 = 3, 6 and 9 JK -1 . The thermostat temperature is kept constant; when the sensor reaches steady state, it increases linearly to +10 K for 150 s and finally remains constant until the end. At the base, a constant power of W1 = 0.3 W is dissipated for 300 s, and when the thermostat temperature increases, it decreases linearly until W1 = 0.1 W, which remains constant

Simulations in the Calibration Base and in the Human Body
The simulations are performed from the systems of equations (Equation (5) or Equation (6)). First, we simulate the operation of the sensor when it is located in its base, which is completed for different heat capacities of the base: C1 = 3, 6 and 9 JK -1 . The thermostat temperature is kept constant; when the sensor reaches steady state, it increases linearly to +10 K for 150 s and finally remains constant until the end. At the base, a constant power of W1 = 0.3 W is dissipated for 300 s, and when the thermostat temperature increases, it decreases linearly until W1 = 0.1 W, which remains constant

Simulations in the Calibration Base and in the Human Body
The simulations are performed from the systems of equations (Equation (5) or Equation (6)). First, we simulate the operation of the sensor when it is located in its base, which is completed for different heat capacities of the base: C 1 = 3, 6 and 9 JK −1 . The thermostat temperature is kept constant; when the sensor reaches steady state, it increases linearly to +10 K for 150 s and finally remains constant until the end. At the base, a constant power of W 1 = 0.3 W is dissipated for 300 s, and when the thermostat temperature increases, it decreases linearly until W 1 = 0.1 W, which remains constant until the end of the simulation. Figure 7 shows the thermostat temperature (T 2 ), the calorimetric signals (y), the power dissipated in the base (W 1 ) and the power dissipated in the thermostat (W 2 ). Additionally, a 0.5 K linear increase in ambient temperature (∆T 0 ) has also been simulated. As expected, a higher value of the heat capacity C 1 implies a slower transient response of the calorimetric signal. until the end of the simulation. Figure 7 shows the thermostat temperature (T2), the calorimetric signals (y), the power dissipated in the base (W1) and the power dissipated in the thermostat (W2). Additionally, a 0.5 K linear increase in ambient temperature (T0) has also been simulated. As expected, a higher value of the heat capacity C1 implies a slower transient response of the calorimetric signal.  Figure 8 shows the simulation of a measurement on the surface of the human body. The simulation takes into account that the sensor is initially at its base (C1 = 3 JK -1 ) and then it is applied to the skin, which has the same or different C1 heat capacity. When placing the sensor on the skin, three phenomena appear that make this measurement different from a calibration measurement at its base: 1) There is an instantaneous contact between the sensor surface and the skin surface that are at different temperatures. This produces a peak in the calorimetric signal, caused by an instantaneous power that is transmitted from the highest temperature surface to the lowest temperature. This instantaneous power that is transmitted from the skin to the sensor is represented with an exponential function of the form where A0 is the power transmitted in steady state and A1 is the amplitude of the decreasing exponential due to instantaneous contact between the sensor surfaces and the skin. In the simulated case represented in Figure 8, we assume a power of A0 = 0.3 W for the initial temperature of the thermostat and a power of A0 = 0.1 W for the final temperature. We assume an amplitude of A1 = 2 W and a time constant of 1 = 9 s.
2) The temperature surrounding the sensor has changed and is not the same temperature surrounding the sensor when it is in the base. In general, the temperature in the neighbourhood of the skin is higher. This temperature difference is represented by T0 and responds to the expression given by (Equation (4)). In the simulated case A = 2.5 K, B = 1 K Figure 7. Simulations for the case of calibrations on a base that has different heat capacities. Curve baselines have been corrected. Figure 8 shows the simulation of a measurement on the surface of the human body. The simulation takes into account that the sensor is initially at its base (C 1 = 3 JK −1 ) and then it is applied to the skin, which has the same or different C 1 heat capacity. When placing the sensor on the skin, three phenomena appear that make this measurement different from a calibration measurement at its base: (1) There is an instantaneous contact between the sensor surface and the skin surface that are at different temperatures. This produces a peak in the calorimetric signal, caused by an instantaneous power that is transmitted from the highest temperature surface to the lowest temperature. This instantaneous power that is transmitted from the skin to the sensor is represented with an exponential function of the form where A 0 is the power transmitted in steady state and A 1 is the amplitude of the decreasing exponential due to instantaneous contact between the sensor surfaces and the skin. In the simulated case represented in Figure 8, we assume a power of A 0 = 0.3 W for the initial temperature of the thermostat and a power of A 0 = 0.1 W for the final temperature. We assume an amplitude of A 1 = 2 W and a time constant of τ 1 = 9 s. (2) The temperature surrounding the sensor has changed and is not the same temperature surrounding the sensor when it is in the base. In general, the temperature in the neighbourhood of the skin is higher. This temperature difference is represented by ∆T 0 and responds to the expression given by (Equation (4)). In the simulated case A = 2.5 K, B = 1 K and a time constant of 9 s are considered. (3) The skin has a heat capacity typical of the area where it is being measured and therefore the transient response will depend on that heat capacity. Figure 8 shows these differences in the calorimetric signal. Three heat capacities for the skin have been used in the simulation: 3, 6 and 9 JK −1 .
Sensors 2020, 19, x FOR PEER REVIEW 11 of 19 calorimetric signal. Three heat capacities for the skin have been used in the simulation: 3, 6 and 9 JK -1 .

Method for Determining Heat Flux and Heat Capacity.
This section describes the calculation method to determine the heat capacity and heat flow on a localized surface of the human body. The effectiveness of this method is verified from the curves simulated in the previous section. The following hypotheses are considered in the simulation of measurements on human skin: 1) When the temperature of the thermostat is constant (T2initial) and the sensor is applied to the skin, the heat flow W1 that passes through the sensor obeys (Equation (9)). Contrastingly, when there is a linear variation in the temperature of the thermostat (from T2initial to T2end), we assume that the power W1 decreases linearly to a final value, which remains constant while the temperature of the thermostat remains constant at its final value (T2end). Thus, the heat flux can be described with (Equation (10)

W t for t t W t A A t t for t t t W t W t A t t t t for t t t W t W t for t t t
In this equation, t1 is the instant in which the sensor is applied to the skin, t2 is the instant in which the linear increase in the temperature of the thermostat begins, t3 is the instant in which the aforementioned linear variation ends and starts to keep the temperature constant, and tend is the final instant of the measurement.

Method for Determining Heat Flux and Heat Capacity
This section describes the calculation method to determine the heat capacity and heat flow on a localized surface of the human body. The effectiveness of this method is verified from the curves simulated in the previous section. The following hypotheses are considered in the simulation of measurements on human skin: (1) When the temperature of the thermostat is constant (T 2initial ) and the sensor is applied to the skin, the heat flow W 1 that passes through the sensor obeys (Equation (9)). Contrastingly, when there is a linear variation in the temperature of the thermostat (from T 2initial to T 2end ), we assume that the power W 1 decreases linearly to a final value, which remains constant while the temperature of the thermostat remains constant at its final value (T 2end ). Thus, the heat flux can be described with (Equation (10)): Sensors 2020, 20, 3431 12 of 20 In this equation, t 1 is the instant in which the sensor is applied to the skin, t 2 is the instant in which the linear increase in the temperature of the thermostat begins, t 3 is the instant in which the aforementioned linear variation ends and starts to keep the temperature constant, and t end is the final instant of the measurement. (2) The difference in ambient temperature ∆T 0 obeys (Equation (4)).
(3) The relationships between all the system variables obey the equations of the model (Equation (5)).
The model parameters have been determined in the calibration (Table 1) except for the value of C 1 which depends on the place of measurement.
The chosen calculation method is similar to that used in model identification and discussed in Section 2.3 of this work. The method consists of determining eight parameters: the first four (A 0 , A 1 , τ 1 , ∆A 0 ) allow us to reconstruct the power W 1 (t). The next three (A, B, τ T ) allow us to build the function ∆T 0 , and the last one (C 1 ) completes the model described by (Equation (5)).
The procedure begins with initial values of the first seven parameters with which the curves W 1 (t) (Equation (9)) and ∆T 0 (t) (Equation (4)) are constructed. From these temporal functions, the power W 2 (t) (known) and the initial value of the eighth parameter (C 1 ), the temperature and calorimetric curve of the thermostat are reconstructed using the system of equations (Equation (5)). From these reconstructions, the error criterion given by Equation (7)) is determined. Using the Nelder-Mead simplex algorithm [33][34][35], the new parameter values are determined. Figure 9 shows a diagram of the calculation procedure. τ1, ΔA0) allow us to reconstruct the power W1(t). The next three (A, B, τT) allow us to build the function ΔT0, and the last one (C1) completes the model described by (Equation (5)).
The procedure begins with initial values of the first seven parameters with which the curves W1(t) (Equation (9)) and ΔT0(t) (Equation (4)) are constructed. From these temporal functions, the power W2(t) (known) and the initial value of the eighth parameter (C1), the temperature and calorimetric curve of the thermostat are reconstructed using the system of equations (Equation (5)). From these reconstructions, the error criterion given by Equation (7)) is determined. Using the Nelder-Mead simplex algorithm [33][34][35], the new parameter values are determined. Figure 9 shows a diagram of the calculation procedure.  Table 2 shows the results obtained for the three simulated measurements represented in Figure   8. The errors in the determination of the heat flow A0 and ΔA0 are less than 0.3%, and in the determination of the heat capacity it is less than 0.4%. The errors (RMSE) in adjusting the calorimetric signal and the thermostat temperature (less than 9 µV and 0.15 mK) are very small in relation to the peak-to-peak value of these curves (120 mV and 10 K). These results on the simulated measurements demonstrate that the proposed method is adequate. Table 2. Results of the calculation on the simulated measurements with C1 = 3, 6 and 9 JK -1 ( Figure  8).   Table 2 shows the results obtained for the three simulated measurements represented in Figure 8. The errors in the determination of the heat flow A 0 and ∆A 0 are less than 0.3%, and in the determination of the heat capacity it is less than 0.4%. The errors (RMSE) in adjusting the calorimetric signal and the thermostat temperature (less than 9 µV and 0.15 mK) are very small in relation to the peak-to-peak value of these curves (120 mV and 10 K). These results on the simulated measurements demonstrate that the proposed method is adequate. Table 2. Results of the calculation on the simulated measurements with C 1 = 3, 6 and 9 JK −1 (Figure 8).

Results
In this section the results corresponding to measurements carried out on two healthy subjects are presented, the characteristic data of both subjects are indicated in Table 3. Both subjects are healthy, although the senior subject suffers from Hashimoto's thyroiditis and is currently receiving treatment.

Measurements in the Junior Subject
Once the suitability of the method in Section 2.5 has been verified, the proposed method is applied to an experimental measurement on human skin. We present a measurement made on the temple of the junior subject. The subject is sitting and resting during the measurement. Four identical consecutive measurements are made. Each of them is of the same type as those performed in the simulation. Figure 10 shows the experimental curves of the calorimetric signal (y), power dissipated in the thermostat (W 2 ) and temperature of the thermostat (T 2 ). Figure 11 shows the placement of the sensor on the subject's temple. The measurements were made in the laboratory with an ambient temperature of 22.6 • C and a relative humidity of 53%.
We apply the calculation method to each measurement. The results obtained are shown in Table 4 where the values of the eight parameters are indicated. Firstly, the parameters (A 0 , A 1 , τ 1 , ∆A 0 ) related to the reconstruction of the heat flow through the sensor as a function of time when applied to human skin. Next, the value of the heat capacity of the area of the human body where the measurement is made is observed, the obtained average value being 5.91 JK −1 . Finally, the parameters (A, B, τ T ) that allow us to reconstruct ∆T 0 , which is the equalization of the temperature around the sensor when it is applied to human skin, are obtained. Adjustment errors are less than 35 µV for a 120 mV peak-to-peak calorimetric signal, and 4 mK for a 10 K peak-to-peak thermostat temperature. Figure 12a shows the heat flow measured by the sensor in the first measurement made on the subject's temple. The initial transient can be observed when applied from the sensor due to the instantaneous contact between the sensor surfaces and the skin. Next, the heat flow reaches a steady state of 288 mW for the thermostat temperature of 26 • C. This power decreases linearly as the thermostat temperature increases, reaching a value of 67 mW for the thermostat temperature of 36 • C. The reconstruction of the calorimetric signal and the thermostat temperature can be seen in Figure 12b, while Figure 12c displays the experimental and calculated curves indicated with the subscripts exp and cal, respectively. The variation ∆T 0 as a function of time is also represented in Figure 12c. temple of the junior subject. The subject is sitting and resting during the measurement. Four identical consecutive measurements are made. Each of them is of the same type as those performed in the simulation. Figure 10 shows the experimental curves of the calorimetric signal (y), power dissipated in the thermostat (W2) and temperature of the thermostat (T2). Figure 11 shows the placement of the sensor on the subject's temple. The measurements were made in the laboratory with an ambient temperature of 22.6 ºC and a relative humidity of 53%.  We apply the calculation method to each measurement. The results obtained are shown in Table 4 where the values of the eight parameters are indicated. Firstly, the parameters (A0, A1, 1, A0) related to the reconstruction of the heat flow through the sensor as a function of time when applied to human skin. Next, the value of the heat capacity of the area of the human body where the measurement is made is observed, the obtained average value being 5.91 JK -1 . Finally, the parameters (A, B, T) that allow us to reconstruct ΔT0, which is the equalization of the temperature around the sensor when it is applied to human skin, are obtained. Adjustment errors are less than 35 µ V for a 120 mV peak-to-peak calorimetric signal, and 4 mK for a 10 K peak-to-peak thermostat temperature. temple of the junior subject. The subject is sitting and resting during the measurement. Four identical consecutive measurements are made. Each of them is of the same type as those performed in the simulation. Figure 10 shows the experimental curves of the calorimetric signal (y), power dissipated in the thermostat (W2) and temperature of the thermostat (T2). Figure 11 shows the placement of the sensor on the subject's temple. The measurements were made in the laboratory with an ambient temperature of 22.6 °C and a relative humidity of 53%.  We apply the calculation method to each measurement. The results obtained are shown in Table 4 where the values of the eight parameters are indicated. Firstly, the parameters (A0, A1, τ1, ΔA0) related to the reconstruction of the heat flow through the sensor as a function of time when applied to human skin. Next, the value of the heat capacity of the area of the human body where the measurement is made is observed, the obtained average value being 5.91 JK -1 . Finally, the parameters (A, B, τT) that allow us to reconstruct ΔT0, which is the equalization of the temperature around the sensor when it is applied to human skin, are obtained. Adjustment errors are less than 35 µV for a 120 mV peak-to-peak calorimetric signal, and 4 mK for a 10 K peak-to-peak thermostat temperature. Figure 11. Sensor application to the skin. Table 4. Results in four consecutive measurements performed on the temple (junior subject, Figure 9).

Model (Equation (5))
∆T 0 (Equation (4))  The power dissipated by the subject depends on the physical situation of the subject, the ambient, the temperature and the humidity, among other factors. We consider that a rigorous thermal characterization of the skin requires, in addition to knowledge of heat flow, data on the heat capacity and thermal conductance of the skin. The heat capacity and heat flux are determined directly by the calculation method. Thermal conductance (or inverse, thermal resistance), can be characterized from the variation of heat flow with the temperature of the thermostat. The inverse of the slope of this line has units of thermal resistance, therefore we define an equivalent thermal resistance of the skin Rskin with the expression given by (Equation (11)). The power dissipated by the subject depends on the physical situation of the subject, the ambient, the temperature and the humidity, among other factors. We consider that a rigorous thermal characterization of the skin requires, in addition to knowledge of heat flow, data on the heat capacity and thermal conductance of the skin. The heat capacity and heat flux are determined directly by the calculation method. Thermal conductance (or inverse, thermal resistance), can be characterized from the variation of heat flow with the temperature of the thermostat. The inverse of the slope of this line has units of thermal resistance, therefore we define an equivalent thermal resistance of the skin R skin with the expression given by (Equation (11)).

Errors (Equation (7))
The thermal resistance of the sensor R sensor is the inverse of the thermal conductivity P 12 . In this case, R sensor = 1/P 12 = 10.41 K/W, ∆T 2 = 10 K, ∆W 1 = ∆A 0 . Substituting the values from Table 4, we obtain R skin = 36.5 ± 1.7 K/W. The uncertainty of the values obtained for the thermal properties of the skin in this series of measurements is 3.5% for heat capacity and 4.7% for equivalent thermal resistance.

Measurements in Senior Subject
In this section, eight measurements performed on different days in the temple of the healthy 62 year-old male are presented. Table 5 shows the results obtained. Figure 13 shows the variation of the heat flow of each day as a function of time and the temperature of the thermostat. The heat flow is very similar on both days, although on the first day it is slightly higher. This could be because the ambient temperature is lower, which produces greater dissipation. As for the thermal properties of the skin in the measured area, the average heat capacity obtained on the first day is 5.7 J/K, on the second day it is 5.9 J/K, and the mean value is C 1 = 5.8 ± 0.2 J/K. The average equivalent thermal resistance obtained on the first day is 31 K/W, on the second day it is 29 K/W, and the average value between the two days is 30 ± 2 K/W. The uncertainty in determining the thermal properties of the skin is 3.4% for heat capacity and 6.7% for thermal resistance. Table 5. Results of measurements made on the temple in the senior subject on two different days (day 1: T room = 23.3 • C, day 2: T room = 23.5 • C).

Model (Equation (5))
∆T 0 (Equation (4)) Errors (Equation (7))   Table 6 shows results of two series of measurements performed on the right hand on two different days. In this case there, is a clear difference in the series of measurements from the first day given that the subject had an unusually low dissipation. The temperature around the sensor in the first day is Troom + A = 21.3 + 0.61 = 21.91 ℃, however, in the second day it is 23.68 ℃. These values clearly show that the surface temperature of the skin is lower in the first case. On the other hand, we can compare the area A1.1 which represents the over-energy of the transitory phenomenon. In the measurements of the first day, said energy is negative or very small compared to that of the second day. This is because the surface temperature of the skin before applying the sensor to the skin is lower on the first day. We also compared the heat capacities and equivalent thermal resistances obtained: on the first day C1 = 5.05  0.20 JK -1 , Rskin = 38.9  0.7 KW -1 , and on the second day C1 = 5.43  0.49 JK -1 , Rskin = 28.9  1.3 KW -1 . The heat capacities obtained in both measurements are less than in the temple and the values obtained in each day are of the same order of magnitude as those that occurred in the measurements made in the temple. Thus, we can consider an average heat capacity in the hand of 5.3  0.4 J/K with an uncertainty of 8%. However, the equivalent thermal resistance obtained in these two measurements is clearly different, so that on the first measured day, with less heat dissipation, the equivalent thermal resistance is 35% greater than that obtained on the second day. With this discussion on the results obtained, we want to indicate that while the heat capacity in an area of the skin keeps its value in a range of  8%, the thermal resistance can change greatly depending on the physical situation of the subject.  Table 6 shows results of two series of measurements performed on the right hand on two different days. In this case there, is a clear difference in the series of measurements from the first day given that the subject had an unusually low dissipation. The temperature around the sensor in the first day is T room + A = 21.3 + 0.61 = 21.91 • C, however, in the second day it is 23.68 • C. These values clearly show that the surface temperature of the skin is lower in the first case. On the other hand, we can compare the area A 1 τ 1 which represents the over-energy of the transitory phenomenon. In the measurements of the first day, said energy is negative or very small compared to that of the second day. This is because the surface temperature of the skin before applying the sensor to the skin is lower on the first day. We also compared the heat capacities and equivalent thermal resistances obtained: on the first day C 1 = 5.05 ± 0.20 JK −1 , R skin = 38.9 ± 0.7 KW −1 , and on the second day C 1 = 5.43 ± 0.49 JK −1 , R skin = 28.9 ± 1.3 KW −1 . The heat capacities obtained in both measurements are less than in the temple and the values obtained in each day are of the same order of magnitude as those that occurred in the measurements made in the temple. Thus, we can consider an average heat capacity in the hand of 5.3 ± 0.4 J/K with an uncertainty of 8%. However, the equivalent thermal resistance obtained in these two measurements is clearly different, so that on the first measured day, with less heat dissipation, the equivalent thermal resistance is 35% greater than that obtained on the second day. With this discussion on the results obtained, we want to indicate that while the heat capacity in an area of the skin keeps its value in a range of ±8%, the thermal resistance can change greatly depending on the physical situation of the subject. Table 6. Results of measurements made on the right hand in the senior subject in two different days (day 1: T room = 21.3 • C, day 2: T room = 21.8 • C).

Discussion
The magnitude of equivalent thermal resistance and heat capacity exposed in this work is global and of non-specific values, and is associated with the measured 2 × 2 cm 2 skin area. The coefficient of variation is 3.5% for heat capacity and 6% for equivalent thermal resistance. These indicators are similar to values obtained in other works [10,15,16]. To compare our values with those in the literature, it is necessary to consider the thermal penetration depth. In our case, this is 3-4 mm depending on the measured area. From this condition, the specific heat capacities and thermal conductivities obtained from our experimental results are consistent with the literature. For example, for a skin volume of 2 × 2 × 0.4 cm 3 , considering an average specific heat capacity of 3.40 Jg −1 K −1 and an average density of 1.15 g/cm 3 [36,37], we obtain an absolute heat capacity of 6.26 JK −1 . For an average thermal conductivity of 0.30 Wm −1 K −1 [36,37], the thermal resistance would be 33.3 KW −1 . These values are similar to our results.
The results obtained in this work have been obtained from experimental measurements made with a previously calibrated instrument. A more precise interpretation of the experimental results requires a better study of the skin's thermal model. An increase in thermal penetration depth implies an increase in the heat capacity and thermal resistance of the zone. Other variables that affect the thermal properties are the physical state of the subject and the presence of lesions in the measurement zone. The physical state of the subject significantly affects the equivalent thermal resistance and, to a lesser extent, the heat capacity.

Conclusions
In this work a calorimetric sensor has been used to measure, by applying it to the skin, the heat flow dissipated by the human body. A non-invasive method is proposed to determine the heat capacity and equivalent thermal resistance of the measured skin area. The modelling of the sensor is effective since it is capable of relating the signals measured by the sensor with the heat flow that passes through the sensor. A sensor calibration has been carried out to determine the invariant parameters of the model. A calculation method that allows one to reliably determine the heat flow and the thermal properties of the skin in the measured area has been proposed. We have verified that the heat capacity of the skin is a property whose value remains at a value that can vary by 8%. However, the equivalent thermal resistance of the skin has a greater variability due to its great dependence on the physical situation of the subject.
With this work, we show the ability of the sensor to determine in vivo the heat flow, the equivalent thermal resistance and the heat capacity of the skin using a non-invasive technique. These physical magnitudes show a great variability and depend on physiological and environmental factors. While this fact could be interpreted as a problem of uncertainty, it is an opportunity to advance the study of human physiology with new data provided by the sensor utilised in this work. The repeatability experiments made in each series show that the variability of the results responds more to changes in the subject and the environment than to changes in the device.
Manually clamping of the sensor on the skin is an inconvenience as it requires specialized personal. Currently, in order to determine the absolute values of the heat flow, it is necessary to transfer the sensor from its base to the skin, and this is carried out manually. However, for the determination of equivalent thermal resistance and heat capacity, it is not necessary for the sensor to be previously on the base; instead, it can be permanently placed on the skin by means of an adapted clamp. We are currently working on this possibility by thermally exciting the skin with sinusoidal thermal dissipations.