Comprehensive impedance spectroscopy equivalent circuit of a thermoelectric device which includes the internal thermal contact resistances

Thermoelectric devices are widely used as solid-state refrigerators and have potential energy generation applications. Their characterization is key to develop more efficient devices and monitor their performance. Electrical impedance spectroscopy has been proved to be a useful method for the characterization of thermoelectric modules. However, deviations from current impedance models still exist in experimental results, especially in the high frequency part of the impedance spectrum, which limits its use. Here, we present a new comprehensive impedance model (equivalent circuit) which covers all the key phenomena that affects the module performance, and it is able to explain the observed deviations. The new equivalent circuit includes, as new additions, the thermal influence of the metallic strips (electrodes), combined with the thermal contact resistance between the metallic strips and the outer ceramic layer. Moreover, a new more accurate spreading-constriction impedance element, which considers the variation of the heat flow in the radial direction at the outer ceramic surfaces, is also developed. The comprehensive equivalent circuit was used to perform fittings to impedance spectroscopy measurements of modules fabricated by different manufacturers. From the fittings, it was possible to identify, among other key properties, the internal thermal contact resistances, whose direct determination is very challenging. Thermal contact resistivities at the metallic strips/thermoelectric elements interface in the range 2.20 × 10-6-1.26 × 10-5 m2KW− 1 were found. An excellent thermal contact was identified at the metallic strips/ceramic layers. This opens up the possibility of using impedance spectroscopy as a powerful tool to evaluate, monitor, and identify issues in thermoelectric devices.


Introduction
Thermoelectric (TE) devices can directly convert heat into electricity and vice versa. They are widely used as solid-state refrigerators and have potential energy generation applications, such as spacecrafts, exhaust gases from industries and automobiles, body heat, solar generators, etc. [1,2]. To evaluate the energy conversion efficiency of TE materials, the dimensionless figure-of-merit zT = S 2 σT/λ is typically employed, where S is the Seebeck coefficient, σ the electrical conductivity, λ the thermal conductivity, and T the absolute temperature. While at a material level different measuring techniques are well established [3], it does not exist an international standard method to characterize TE devices [4]. Solving this widely recognized problem is critical, as device-level characterization is needed to validate the implementation of novel TE materials [5], guide system design [6], perform device quality control, and monitor device performance and its possible degradation [4]. It is therefore crucial for the development of the TE technology to achieve standardized, accurate and easy-to-implement characterization methods.
Current measurement methods have several disadvantages which limit their use: they usually require intricate and expensive experimental setups, the use of multiple measuring systems, and are known to potentially have large uncertainties [4,7]. In recent years, electrical impedance spectroscopy has been successfully used to measure the dimensionless figure of merit zT of TE modules [8][9][10][11]. In addition, the possibility of performing a complete characterization of TE devices (determination of the internal ohmic resistance, the average Seebeck coefficient and thermal conductivity of the TE legs, and zT) under suspended conditions if the thermal conductivity of the ceramics is given, has been demonstrated [12][13][14][15], without the need of multiple measuring systems. To perform the complete characterization, an ideal equivalent circuit was proposed by us [12], which was obtained by solving the heat equation in the frequency domain. Obtaining the system equivalent circuit is important, since once known, characterization can be performed by fitting the experimental data to the equivalent circuit. In that study [12], the TE device was thermally modelled by considering multiple TE legs in contact with outer ceramic layers of same cross-sectional area, which led to an equivalent circuit formed by a resistor in series with the parallel combination of constant-temperature and adiabatic Warburg elements.
The ideal model was improved by adding spreading-constriction effects [13,16], which take into account the differences in the crosssectional area between the TE legs and the outer ceramic layers. Another improvement of the ideal model was obtained by adding the convection effect at the outer ceramic surfaces [17]. Three new elements in the equivalent circuit were obtained due to this effect. A more comprehensive equivalent circuit included the radiation effect [18] in conjunction with the spreading-constriction and the convection phenomena [13]. Also, it has been included recently the effect of the metallic strips (electrodes), assuming they behave like a capacitor, and the presence of a thermal contact resistance between the thermoelements and the electrodes [19]. Apart from these models under suspended conditions, we have recently reported an electrical impedance equivalent circuit of a module attached to heat sinks at both sides, which takes into account the effect of a thermal contact resistance between the outer ceramic surfaces and the heat sinks [20].
It is worth mentioning that the spreading-constriction impedance was previously reported [13] considering that the heat flux at the outer ceramic surfaces is uniform, which is only true for TE modules suspended in vacuum (without convection) and at room conditions (assuming radiation negligible). An alternative spreading-constriction impedance expression was obtained [16] considering that the temperature at the outer ceramic surfaces is uniform, but only a perfect contact with an ideal heat sink could achieve such condition, and this expression also fails to cover intermediate cases.
In this work, we present experimental electrical impedance measurements performed at different TE modules that show features, specially at high frequencies, which cannot be explained by current impedance models. This limits the possibility of widely using this highlybeneficial method for TE module assessment. To avoid this, and to explain the different features experimentally observed, we have developed a new more comprehensive equivalent circuit that includes all the key phenomena affecting the TE module performance. Namely, the effect of the thermal contact resistance between the TE legs and the metallic strips that connect them, the own electrodes contribution, and the thermal contact resistance between the electrodes and the outer ceramic plates. Moreover, a new spreading-constriction impedance element, which considers the variation of the heat flux in the radial direction at the outer ceramic surfaces, is obtained to cover non-ideal cases. This new equivalent circuit, which also includes convection and radiation effects previously developed, and the inductive phenomena existing at the highest frequencies [21], allows the analysis of the different characteristic features experimentally observed, which was not possible before with the previous models. Furthermore, the new equivalent circuit was also used to perform fittings to TE modules from different manufacturers to determine their internal thermal contact resistance. The fact that impedance spectroscopy can identify and quantify all these phenomena, in addition to its capability to accurately measure the module zT from a single and simple measurement, opens up the possibility of using this method as a standard quality control tool, able to identify and monitor in detail issues in TE modules in energy applications.

Theoretical model
In order to obtain an electrical impedance model which includes all the key phenomena that can be relevant in standard TE devices, the model shown in Fig. 1 was considered, which consists of a cylindrical TE leg of cross-sectional area A and length L, contacted by two metallic strips (usually copper) with a slightly larger cross-sectional area A/η M (being η M the ratio between the area of all the TE legs and the area of all the metallic strips) and length L M , and two external ceramic layers of cross-sectional area A/η (being η the filling factor of the TE module, i.e. the ratio between the area of all the TE legs and the ceramic area, typically around 0.3), and length L C . The use of cylindrical legs has been previously adopted [13,16] and has also shown to have no significant differences with respect to a prismatic geometry [22], and it simplifies the thermal spreading-constriction analysis.
At both sides of the metallic strips, the possibility of having a thermal contact resistance is considered, including r TC1 and r TC2 as thermal contact resistivities between the TE legs and the metallic strips, and between the metallic strips and the ceramic layers, respectively. In addition, radiation/convection effects are included around the TE legs (h), around the metallic strips (h 1 ), at the inner ceramic surfaces (h 2 ), and at the outer ceramic surfaces (h 3 ), as shown in Fig. 1 [13,17,18], being approximately, where σ B is the Stefan-Boltzmann constant, and ε TE , ε M , and ε C the average equivalent emissivity of the TE legs, metallic strips, and external ceramic materials, respectively. The radiation expression assumes small temperature increase from ambient temperature, which is the case of impedance spectroscopy measurements, and has previously been employed in similar derivations [18]. Moreover, h ic and h ec are the average internal and external convection heat transfer coefficients, respectively, and are assumed to be constant. It is worth noting that the internal (-L M ≤ x ≤ L + L M ) convection may be negligible in some cases (e.g. large amount of TE legs in the TE module or use of sealants), and the radiation effect depends strongly on temperature.
Furthermore, a TE module sandwiched between heat exchangers that do not change their temperature (ideal heat sinks) with a thermal contact resistivity r TC can also be considered at the outer ceramic surfaces instead of convection/radiation, being in this case h 3 = 1/r TC [20]. Joule effect is neglected since the ac amplitude used in electrical impedance spectroscopy measurements is small and TE materials usually have relatively high electrical conductivity. It should also be noted that all the TE properties are considered temperature independent, and that the final response must be multiplied by the number of legs, 2 N (being N the number of TE couples), since it is assumed that all the TE legs are identical.
To account for the spreading-constriction impedance introduced at x = -L M , and at x = L + L M by the change of the cross-sectional areas, we have developed a new expression considering that the heat released/ absorbed at the outer external ceramic surfaces varies radially. It should be noticed that the spreading-constriction at x = 0 and at x = L can be neglected due to the typically high thermal conductivity of the metallic strips [16]. The new expression was obtained following a similar procedure used by Casalegno et al. [16], with the exception that no constant temperature but convection and radiation at the outer ceramic surfaces was assumed, where θ 3 is the temperature with respect to the initial temperature T initial in the frequency domain [θ=L (T-T initial )] at the outer ceramic surfaces, which may change radially. This change leads to a slightly larger expression for the spreading-constriction impedance (see its derivation in Annex I), with λ C being the thermal conductivity of the ceramic layer, J 0 and J 1 the first kind Bessel functions of order zero and one, respectively, r M and r C the equivalent radii of the metallic strip and ceramic layer, respectively, δ n is the nth zero of J 1 , and γ n is the value for each δ n that verifies, where j=(-1) 0.5 is the imaginary number, ω the angular frequency (ω = 2πf, being f the frequency), and α C the thermal diffusivity of the external ceramic material. It is worth highlighting that Eq. (6) agrees with the solution obtained for constant heat flux at the outer external ceramic surfaces [13] when the convection/radiation effect is negligible (h 3 → 0). It also agrees with the solution obtained for constant temperature condition on the same surfaces [16] when a perfect contact with a heat exchanger is considered (h 3 →∞). Finally, this expression agrees with the solution obtained by Yovanovich et al. [23] for steady state conditions (ω → 0).
The electrical impedance Z = V/I of a TE module is given by, where V(0) and V(L) are the voltages at x = 0 and x = L, respectively, I 0 is the electrical current crossing the device at x = 0, R Ω is the total ohmic resistance of the TE module, which includes the contribution of all the TE legs, the metallic strips, the wires, and the electrical contact resistances, S is the average Seebeck coefficient of the thermoelements, and T(0) and T(L) are the temperatures at x = 0 and x = L, respectively. It is worth noting that in Eq. (8) the temperature difference between the thermoelement sides can be determined from the temperature at x = 0, due to the anti-symmetry of the system at x = L/2 (see Fig. 1).
To determine T(0) the two-dimensional heat equation in the frequency domain for the three layers (thermoelements, electrodes and ceramics) should be solved, which in cylindrical coordinates takes the form, where r and × are the radial and axial axis (shown in Fig. 1), respectively, and α i the average thermal diffusivity of each material: TE legs (i = TE), metallic strips (i = M), and ceramic layers (i = C). However, averaging the temperature distribution over the equivalent circular surface (θ) of radius r k (being k = TE for the TE legs, k = M for the metallic layers, and k = C for the ceramic layers), the two-dimensional heat equation given in Eq. (9) can be approximated to a one-dimensional heat equation [24], with λ i being the average thermal conductivity of each material, TE legs (i = TE), metallic strips (i = M), and ceramic layers (i = C), and h i the convection/radiation effects around the TE legs (h), and around the metallic strips (h 1 ). Notice that the convection/radiation from the edges of the ceramics can be neglected and hence the term 2h i /(λ i r k ) is not included for this layer (see Fig. 1).
We can now make use of the definition of the total thermal admittance y 0 = ϕ 0 /θ 0 , being ϕ 0 the heat power at x = 0 and θ 0 = θ(x = 0) at the TE leg side. The total thermal admittance, is the series summation of the thermal admittance towards the TE element, y TE , and towards the metallic strip, y e [9], Since the heat power at x = 0 is the Peltier heat [due to the small ac current amplitude applied, it is approximately ϕ 0 =-|S|T initial i 0 , being i 0 =L (I 0 )], Eq. (8) can be rewritten in the frequency domain as, For the determination of y TE and y e , the thermal quadrupole method was used [24]. In this method, the different equations are displayed in different matrices, where each matrix represents one layer (TE material, metal strip, or ceramic) or one thermal restriction (thermal contact resistances, spreading-constriction impedance, and convection/radiation effects). The matrix that defines the TE elements admittance (1/y TE = θ 0, TE /ϕ 0,TE ) is only considered until the half length of the TE layer, since a boundary condition at this position (x = L/2) exists due to the anti- where ω TE is the characteristic angular frequency [ω TE = α TE /(L/2) 2 ], being α TE the average thermal diffusivity of the TE legs, λ TE the average thermal conductivity of the TE legs, and ϕ L/2 the heat flow in the frequency domain at x = L/2.
The matrix that defines the heat flow going towards the metallic strip (1/y e = θ 0,e /ϕ 0,e ) is given by seven matrices, which must be written in a proper order. The first matrix corresponds to the thermal contact resistance between the TE layer and the metallic strip. The second matrix corresponds to the metallic strip itself. The third one to the thermal contact resistance between the metallic strip and the ceramic layer. The fourth matrix is the spreading-constriction impedance of the heat entering the ceramic (which has a larger area than the metallic strip). The fifth one denotes the heat lost in the inner part of the ceramic due to convection/radiation. It should be noticed that in this matrix the area of ceramic considered to be exposed to inner convection is slightly higher than it should, since it also includes the exposed part of the metal (r Mr TE ), which is adopted for simplicity. The sixth one relates to the ceramic layer. Finally, the seventh matrix corresponds to the boundary condition at x = -L M -L C , which describes the convection/radiation at the outer ceramic surfaces (ϕ 3 = θ 3

h 3 A/η),
where ω M = α M /L M 2 (being α M the thermal diffusivity of the metallic layers), ω C = α C /L C 2 , λ M and λ C are the characteristic angular frequencies and thermal conductivities of the metallic strips and ceramic layers, respectively.
Once the two admittances are defined, they are introduced in Eq. (13) using Eq. (12), to obtain the electrical impedance function. A [ parasitic inductance impedance (jωL p ) is also added to this impedance function, which accounts for the inductive response at the highest frequencies, being L p the inductance [21]. After some algebraic steps, the equivalent circuit shown in Fig. 2a is obtained and takes the form, where R Ω is the total ohmic resistance of the TE device, and Z TOT1 and Z TOT2 are defined as, The elements in Eqs. (16), (17), and (18) are defined by, This equivalent circuit contains the total ohmic resistance (which includes contributions from the wires, the metallic strips, the TE legs, and the electrical contacts) of the TE device (R Ω ), the constant-temperature Warburg element (Z WCT ) due to the TE legs, and the adiabatic Warburg element (Z Wa ) due to the external ceramics, which are the three elements of the ideal model [12]. It also contains the resistance R h3 = 4NS 2 T initial η/(h 3 A), the constant-temperature Warburg impedance Z WCT,C , and the capacitor C h3 = R h3 /(R C 2 ω C ). These 3 elements appear when heat is exchanged at the outer ceramic surfaces of the TE module (either by convection/radiation in suspended modules, or by the effect of heat sinks when they are contacted to the module) [13,17,20]. (1-η)], due to the heat removal on the inner ceramic surfaces, and the spreading-constriction impedance Z S/C , due to the effect of the area variation between the metallic strips and the external ceramics, are also present, as previously reported [19]. However, it should be noted that the position of R h2 identified in the equivalent circuit in our analysis (in parallel with the combination of the four elements of the ceramics (Z Wa , R h3 , Z WCT,C , and C h3 ) and in series with Z S/C , see Fig. 2) is different from that shown in [19], where it is in parallel with Z S/C . This is due to a different position adopted here for the matrix of the spreading-constriction [between the z s/c and the external ceramic matrices, see Eq. (15)] which is more correct, since the heat losses from the internal surface of the ceramics is not possible if the heat flow is not spread. In any case, it is expected that this modification will not produce large deviations for commercial TE modules in suspended conditions, since the heat removal by the internal convection and radiation effects is low and the thermal conductivity of the ceramics is high.
The resistance R TC1 = 4NS 2 T initial r TC1 /A, due to the presence of a thermal contact resistance between the TE legs and the metallic strips, is also part of the equivalent circuit of Fig. 2. When this thermal contact exists, part of the heat flowing towards the ceramic is blocked. If this thermal contact is very large (R TC1 →∞), the heat only flows towards the TE legs, and only the Z WCT element (in series with L p and R Ω ) is observed. On the other hand, if this thermal contact is negligible, R TC1 = 0, this element can be replaced by a short circuit (graphically it will be represented by a straight line in Fig. 2). The resistance R TC1 and Z Wa,M , which can become a capacitor for metallic strips of high thermal conductivity, as it is the case for copper [25], were previously identified when the effect of the metallic strips was taken into account in the impedance response [19].
In addition to all the aforementioned elements, nine new elements appear: two resistances (R TC2 and R h3,M ), one adiabatic Warburg element (Z Wa,C,M ), two constant-temperature Warburg elements (Z WCT,M and Z WCT,C,M ), a new expression for the spreading-constriction impedance (Z S/C,M ), and three capacitors [ . The resistance R TC2 = 4NS 2 T initial r TC2 η M /A comes from the thermal contact resistance between the metallic strips and the ceramic layers, and it appears twice in the equivalent circuit since it is also present in the denominator of the capacitor C TC2 [see Eq. (30)]. R TC2 = 0 and becomes a short circuit if the thermal contact is negligible, but for a large thermal contact resistivity (R TC2 →∞), it blocks all the heat flow towards the ceramic elements, leaving the contribution to the equivalent circuit of both the metallic strips and the external ceramics with only the adiabatic Warburg element of the metallic strips (Z Wa,M ) in series with R TC1 , since Z CTC2 → 0. A more detailed information about the physical meaning of resistors and capacitors that come from the existence of a thermal contact resistance can be found in ref. [20].
The elements Z Wa,M , and Z WCT,M appear in the equivalent circuit due to the metallic strips, as it was the case in our previously published article when we added the effect of the convection at the outer ceramic surfaces to the ideal model [17]. The element C h2 = λ M ρ M C p,M Aη/ [4NS 2 T initial ηh 2 (1-η)η M 2 ] is a capacitor that it is influenced by two thermal parameters, h 2 and the thermal effusivity of the metallic strips e M =(λ M ρ M C p,M ) 0.5 , Both h 2 and e determine the temperature of the junction at the side of the metallic layer, which is eventually governed by the heat release/accumulation at the interface [20]. If the metallic strips are large, the accumulation of heat in C h2 will be high, but that accumulation will be lower when heat losses through the internal surface of the ceramics increases (higher h 2 ). The elements Z S/C,M , Z WCT,C,M , C h3,M , Z Wa,C,M , and R h3,M appear in the equivalent circuit when the influence of the metallic strips cannot be neglected. It should be noticed that all these elements involve the properties of the metallic strips, and they disappear if the metallic layers are not present (L M = 0).
It should be noted that when the TE modules are measured suspended under vacuum (no convection), and the radiation is negligible, h, h 1 , h 2 , and h 3 become zero, and the simplified equivalent circuit of Fig. 2b is obtained. This simplified circuit is easier to use at the time of performing fittings to extract the parameters of the system.

Results and discussion
In order to evaluate the capabilities of the new comprehensive equivalent circuit developed, the impedance response of three commercial Bi-Te TE modules from different manufacturers were measured:  Table 1. An I ac = 30 mA was used for the three TE modules after its optimization as described in ref. [26]. Amplitude optimization basically consists in identifying the lowest possible current amplitude (I ac ) with no noise in the spectra. The frequency range from 10 mHz to 1 MHz and 50 measuring points (logarithmically distributed in the frequency range) were chosen to ensure a proper number of points in the regions of interest in the spectra, mainly in the high frequency part. All the measurements were performed with the modules suspended under vacuum (<5x10 -4 mbar) at room temperature with a PGSTAT302N potentiostat (Metrohm Autolab B. V.) equipped with a  FRA32M impedance module. Dots in Fig. 3a, Fig. 4a, and Fig. 5a show the experimental impedance spectrum of Module 1, Module 2 and Module 3, respectively. It can be observed that they all differ in their response at high frequency (inset of each figure). Fittings to these experimental results (lines in Fig. 3a,  Fig. 4a, and Fig. 5a) were performed using the equivalent circuit of Fig. 2b by means of our own Matlab code (provided in the Supplementary Information). When using the code, it is possible to choose the frequency range and the parameters to fit. Also, the code can be used to simply perform simulations. In the Supplementary Information, the experimental data of Module 3 is also provided as an example of experimental data. It should be noted that some points at the highest frequencies (4 points for Module 1 and Module 2, and 6 points for Module 3) were not fitted since they do not behave purely inductive.
To obtain the fittings for the different modules using our Matlab code, we followed several steps. First, the specifications of each of the TE modules (see Table 1) were provided. Then, from all the different parameters that can be fitted (L p , R Ω , r TC1 , r TC2 , S, λ TE , α TE , λ C , α C , λ M , and α M ), we selected the key ones to fit, and provided fixed values for the others. When the aim is to fit the experimental data to obtain the values of the internal thermal contact resistivities (r TC1 and r TC2 ), it is recommended to fit L p , R Ω , r TC1 , r TC2 , λ TE , and λ C , and maintain fixed S, α TE , α C , λ M , and α M . It should be taken into account that since many processes overlap at the high frequency region, it will become very difficult to fit the experimental data leaving all the variables free in most of the cases, since the effect of the variation of one parameter can be vanished by the variation of another, and the fittings will converge in unrealistic solu- When the fittings were performed in this way, the value of r TC2 resulted in the order of 10 -7 m 2 KW − 1 , which is quite low, and with an error > 100% for the three TE modules evaluated. This indicate that such thermal contacts are really good and can be neglected. This is not a surprising result, since typically the most problematic junction in TE devices is the one formed by the metallic strips and the thermoelements [4,27], and the thermal contact between the metallic strips and the ceramic layer is more robust and less problematic. After reaching this result, we neglected r TC2 (it became fixed with a zero value) in a second fitting to all the modules. The results of these final fittings can be found in Table 2.
All the values obtained in Table 2 are similar to typical values found in the literature [7,13,28,29]. It should be noted that r TC1 can be directly obtained from the fittings, which is very challenging to be determined by other methods. Moreover, this can be achieved from a single measurement that only takes a few minutes and only involve suspending the module under vacuum and the measurement of its Seebeck coefficient. We would like to remark that the measurement of exceptionally good thermal contacts (contacts with very low thermal contact resistance) and/or different TE module geometries may produce smaller changes in the impedance spectra, which may be more difficult to fit. Furthermore, the values obtained from the fitting are the average properties of the entire TE module, so small local changes might not be sensed, specially of thermal parameters. It should be also noted that the developed equivalent circuit will not be adequate if deviations from the assumptions made in the theoretical model (e.g. significant differences in the TE properties of n and p-type legs, non-uniform leg cross-sectional area, and cylindrical module architectures) exist.
It should be noticed that the r TC1 values are similar to the predicted values obtained using a qualitative method previously developed by us in a PhD thesis [30], although the methodology described in the thesis did not include the inductance contribution, which introduces deviations. Fig. 3b, Fig. 4b, and Fig. 5b show simulations which include the spectrum from the fitting to the experimental data and spectra simulated using the same parameters from the fitting but varying the r TC1 value. Fig. 3c, Fig. 4c, and Fig. 5c compare the fitted spectra with the simulated response without the inductance (L p = 0). Finally, Fig. 3d, Fig. 4d, and Fig. 5d include simulations which include the spectrum from the fitting to the experimental data and spectra simulated using the same parameters than the fitting but varying the r TC2 value.
For Module 1, it can be observed from the inset of Fig. 3b that a

Table 2
Fitting parameters with their associated relative errors (in brackets) obtained for the three TE modules used in this study. The fittings were performed with the Matlab code provided in the Supplementary Information. variation of r TC1 towards higher values leads to smaller slopes and to a more linear trend, which also results in a larger size of this part of the spectra. In addition, as r TC1 is increased, it can be observed from Fig. 3b that the intercept with the real axis occurs at higher real impedance values, due to the overlap with the inductive part, which becomes slightly bent as r TC1 increases from 0 to 2.2x10 -6 m 2 KW − 1 . In Fig. 3c, it is shown the effect of the inductance. Its presence pushes the initial features of the impedance spectrum that appears at the highest frequencies into the Z''/Z' quadrant. When this occurs, a slight tilting in the impedance response is observed before reaching the real axis (Z''=0), instead of having a purely vertical trend. Finally, for Module 1, it can be observed from the inset of Fig. 3d the effect of the presence of r TC2 , which increases the size of this part of the spectra, as occurred for r TC1 as well (see inset of Fig. 3b). Moreover, as r TC2 increases a more curved trend is observed, which is not the case for the experimental spectrum, which shows a more linear trend in this part (inset of Fig. 3a). For this reason, r TC2 was discarded by the fitting (relative errors > 100%) as it was mentioned above.
A similar analysis to that performed for Module 1 can be applied to Module 2. As observed from the inset of Fig. 4b, higher values of r TC1 result in a wider and more horizontal-like response, and the intercept with the real axis is again shifted by the presence of r TC1 . From Fig. 4c, it can also be seen that the inductance pushes the initial features of the impedance spectrum into the Z''/Z' quadrant. Finally for Module 2, it can be observed from the inset of Fig. 4d that the effect of the presence of r TC2 induces a curvature and a larger size in the response of this part of the spectra, which is not the case for the experimental spectrum that exhibits a more horizontal-like trend in this part (see inset of Fig. 4a). This led to the rejection of r TC2 by the fitting. It should be noted from the analysis performed for modules 1 and 2 that even though the length of the metallic strips of these modules is short (see Table 1), the impedance method is able to distinguish between r TC1 and r TC2 . Fig. 5 shows the experimental impedance response and simulations for Module 3. In the experimental response, it is observed from the inset of Fig. 5a a somewhat different pattern compared to the other two modules. An initial curvature is followed by a straight line. As it can be seen in the inset of Fig. 5b, the slope of the straight-line region decreases with r TC1 , at the same time that this part of the spectra becomes larger. From the inset of Fig. 5c, it can be seen that the inductance pushes again the initial features of the impedance spectrum into the Z''/Z' quadrant, and in this case reduces significantly the size of the first curvature. Regarding the influence of r TC2 , it is observed from the inset of Fig. 5d that r TC2 introduces a curvature in the straight-line region, which is not observed experimentally and hence was also rejected by the fitting as mentioned above.

Conclusions
A comprehensive equivalent circuit, which covers all the key phenomena that affects the performance of thermoelectric modules, has been developed. The new equivalent circuit includes, as new additions, the thermal influence of the metallic strips, combined with the thermal contact resistance between the metallic strips and the outer ceramic layer. Moreover, a new more accurate spreading-constriction impedance element, which considers the variation of the heat flow in the radial direction at the outer ceramic surfaces, was also developed. Apart from these new additions, the equivalent circuit includes other key different elements previously reported. The developed equivalent circuit was used to fit experimental measurements of three commercial Bi-Te modules from different suppliers and their internal thermal contact resistances were obtained. To perform the fittings, a new Matlab code has been created and made available (Supplementary Information). The fitting results revealed the presence of thermal contact resistances at the metallic strips/thermoelectric elements interface in all the commercial modules, with thermal contact resistivity values of 2.20 × 10 -6 , 5.29 × 10 -6 and 1.26 × 10 -5 m 2 KW − 1 for the three modules. However, no significant influence of the metallic strips/ceramic layers thermal interface was identified. Moreover, it was found that the inductance of the system plays a crucial role at the highest frequencies, which may produce large deviations in the results if not considered. This study opens up the possibility of using impedance spectroscopy as a powerful tool to detect and measure thermal contact resistances inside thermoelectric devices, which is very challenging for other techniques, and monitor in detail possible issues that could appear during their performance or manufacturing process.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.