Although some studies propose the usage of methane (not natural gas) mixed with hydrogen as an alternative to fossil fuels (Makaryan et al., 2022), we aim to explore the feasibility of a sensing system for mixtures where the higher hydrocarbons and other elements cannot be neglected—if pure methane is mixed with pure hydrogen, a simple physical detection system can be readily proposed. Before proceeding with our analysis, it is instructive to observe the composition of several typical natural gas samples extracted from different sources, shown in Table 1. Notice that some gases are classified as “H” and others as “L”, denoting a high or low percentage of methane, respectively.

One approach to characterising hydrogen-enriched natural gas would be the determination of the molar percentages of its most important components. A rigorous PGC-like, molecule-based approach would require the computation of nine unknowns, which can only be achieved if nine independent equations are available (one of them being given by the definition that the sum of the molar fractions is one). Recognising that the combination of pentane, hexane and higher hydrocarbons sums up less than about 0.1 mol%, their contribution can be either ignored or attributed to a slightly higher amount of butane without causing a significant error. Even so, 7 M percentages remain to be determined, which requires at least six independent physical measurements.

The difficulty in determining the full composition of this complex mixture has led some to search for correlation tendencies (Huber et al., 2018). What is of interest to the application is not the molar percentage of every component, but only that of hydrogen, alongside the overall combustion enthalpy. Hence, if an empirical combination of physical variables can be found, such that the two sought-after variables can be predicted within some accuracy, the computation of the composition can be avoided. Such a finding cannot be guaranteed, as it does not follow from a rigorous model, but it would be of great practical value if it is confirmed. In this article, we test some of the proposed correlation tendencies, finding negative results. The possibility of developing a machine-learning algorithm to predict these two variables based on a certain set of measurable properties, however numerous, remains an open possibility to be explored.

We address the question of a plausible physical characterisation with a different approach, considering hydrogen-enriched natural gas, not as a mixture of gas molecules, but as a mixture of mixtures. This is illustrated in Figure 1. Under this approach, each source of natural gas is characterised by a fixed composition and therefore by a set of effective thermodynamic properties that could serve as its fingerprint. If a certain user draws natural gas from two different sources and hydrogen is added to it, the quantification problem reduces to solving a system of two equations with two unknowns. This assumption is only valid for ideal gases, where molecule-molecule interactions are neglected. Therefore, our approach is limited to pressures less than about 10 bar (Cézar de Almeida et al., 2014) and so not applicable to ‘high-pressure’ pipelines, which are designed for pressures above 16 bar (Télessy et al., 2024). In the case of real-gas behaviour, further variables and more complex equations are required to characterise the mixture.

With a mixture-based approach, a system of equations of the following sort (Equation 1) is solved: x H + x 1 + ⋯ + x n = 1 f i x H , x 1 , … , x n = y i … (1) where x H is the amount fraction of H 2 , x 1 … n the amount fractions of the different ‘kinds’ of natural gas, and y i is a physical property that can be explicitly calculated as a function of the molar percentages. Solving for all the involved molar percentages, not just that of hydrogen, is important in order to estimate the enthalpy of combustion (Equation 2). Δ H C = x H Δ H H + ∑ k = 1 n x k Δ H k (2)

If a sensor system is capable of estimating the enthalpy of combustion, the rate of energy transported through the gas can be then computed according to Equation 3, assuming ideal gas behaviour. (Let v be the flow velocity, A the area of the cross section of the conduit, P the absolute pressure, R the ideal gas constant, and T the absolute temperature.) In other words, a sensor system for billing purposes requires a measurement of flow velocity, absolute pressure, and temperature. E ˙ = n ˙ Δ H C = v A P R T Δ H C (3)

The ease or difficulty of solving this problem lies in the nature of the equations for each physical property, y i . There ought to be an explicit equation under which each property can be directly related to each amount fraction. The equations need not be linear nor polynomial, but any valid approximation as either of these would significantly reduce the complexity.

Let us now consider the available models through which the aforementioned thermodynamic properties (density, speed of sound, thermal conductivity and viscosity) can be calculated for a mixture of gases.

The ideal gas equation yields a simple linear model to calculate the density of a gas mixture (Equation 4) Rather than density itself ( ρ ), it is more convenient to refer to the molar mass ( M ) of the mixture, for it is independent on the absolute pressure and temperature. A measurement of density, therefore, ought to be accompanied by pressure and temperature sensors in order to make the respective compensation. This molar mass is simply the linear combination of the individual molar masses and the respective molar percentages. Notice that this linear combination can be expressed both in terms of the individual molecules and in terms of the different sources of natural gas. ρ R T P = M = ∑ M i x i (4)

The Newton-Laplace equation for the speed of sound ( c ) in an ideal gas reveals its dependence on the overall molar mass of the mixture ( M ), as well as on the overall molar isobaric heat capacity ( C p ) and molar isochoric heat capacity ( C v ), all of which are given by linear combinations of the individual amount fractions (Dixon and Greenwood, 1925). Given the temperature dependence of the speed of sound, it is more convenient to calculate the ratio R T / c 2 , which has units of molar mass and is denoted here as M * (Equation 5). The presence of linear combinations for the calculation of C p , C v and M also makes it possible to express the speed of sound in terms of either the molar percentages of the individual molecules or of the different ‘mixtures’ in the mixture (each with effective values of C p , C v and M ). R T c 2 = M * = ∑ C v i x i ∑ M i x i ∑ C p i x i (5)

Closer examination of Equation 5 leads to the observation that a quadratic polynomial of the individual amount fractions is always formed. For the case of hydrogen mixed with two different compositions of natural gas ( x 1 and x 2 , with x 2 = 1 − x H − x 1 ), Equation 5 can be rewritten as the quadratic polynomial in Equation 6, whereby each coefficient ( ζ i j ) originates from a multiplication between a heat capacity and a molar mass. ζ 20 x H 2 + ζ 11 x H x 1 + ζ 02 x 1 2 + ζ 10 x H + ζ 01 x 1 + ζ 00 = 0 (6)

The calculation of the thermal conductivity and viscosity of a gas mixture—even of ideal gases—is not a straightforward task. Mason and Saxena (1958) developed an approximate formula for the thermal conductivity based on a rigorous analysis of the kinetic theory of gases. A similar work was performed by Wilke (1950) for the calculation of viscosity. Both equations are very similar in form, showing how the pairwise interactions between components affect the overall behaviour of the mixture.

The thermal conductivity of a mixture ( κ m i x ) is thus calculated according to ³ Equation 7: κ m i x = ∑ i = 1 n κ i 1 + 1 x i ∑ j = 1 , j ≠ i j = n x j G i j G i j = 1.065 2 2 1 + μ i μ j 1 2 M j M i 1 4 2 1 + M i M j 1 2 (7)

And the viscosity of a mixture ( μ m i x ) is calculated according to Equation 8: μ m i x = ∑ i = i n μ i 1 + 1 x i ∑ j = 1 , j ≠ i j = n x j ϕ i j ϕ i j = 1 + μ i μ j 1 2 M j M i 1 4 2 4 2 1 + M i M j 1 2 (8)

Both models assume that the individual thermal conductivities and viscosities of the components are known, and that they are all referred to the same temperature. A model to describe the dependence of these properties on the state variables is further required to make the corresponding compensation. A first approximation for ideal gases can be drawn from kinetic theory, whereby the thermal conductivity varies with the square root of temperature and is independent of pressure (Mason, 2023), similarly to the speed of sound.

The reader may notice that the form of Equations 7, 8 gives little opportunity for a closed-form solution of a system of equations. Such a system can be solved numerically with the aid of search algorithms, but it could be impractical for an implementation in a portable electronic device. Here we apply and test a second-order Taylor approximation, knowing that the definition of the problem limits the concentration of hydrogen to relatively small values and that natural gas is mostly composed of methane. We test this polynomial approximation against rigorously calculated properties calculated according to Equations 7, 8. This second-order approximation has the simpler form in Equation 9 for the case of a mixture of hydrogen plus three different sorts of natural gas ( x 1 , x 2 , x 3 with x 3 = 1 − x H − x 1 − x 2 ). κ m i x ≈ κ 200 ⁢ x H 2 + κ 020 ⁢ x 1 2 + κ 002 ⁢ x 2 2 + κ 110 ⁢ x H ⁢ x 1 + κ 101 ⁢ x H ⁢ x 2 + κ 011 ⁢ x 1 ⁢ x 2 + κ 100 ⁢ x H + κ 010 ⁢ x 1 + κ 001 ⁢ x 2 + κ 000 (9)

In order to verify the accuracy of the proposed estimation methods, a set of randomly generated—yet feasible—natural gas compositions were calculated. Table 2 shows the tabulated properties of the relevant molecules at 20°C and atmospheric pressure, extracted from the NIST database (Lemmon et al., 2023). With this information, the thermodynamic properties of several random mixtures of natural gas and hydrogen can be computed, using the models in Equations 3, 4, 6, 7. We ‘synthesised’ thus one hundred random binary combinations of the seven natural gas sources in Table 1 and then calculated their properties upon mixture with hydrogen from 0 mol% to 30 mol%. If the proposed estimators are able to compute the combustion enthalpy and the hydrogen content of such a wide spectrum of natural gas compositions, such a system can certainly be calibrated for a more regulated application case.

In this study, we compare four kinds of estimators with increasing complexity. The simplest estimator is based on a measurement of the thermal conductivity alone and is calibrated against one source of natural gas (“NG type 1” in Figure 2). It naïvely assumes that there is no variation in the composition of natural gas, and any change in the thermal conductivity is attributed to a change in the hydrogen content alone. The combustion enthalpy is then simply calculated as the linear combination of the enthalpy of hydrogen and that of the calibration gas with the corresponding proportions. The following two estimators are based on two variables ⁴ —speed of sound together with thermal conductivity, and density together with thermal conductivity. A two-variable estimator requires a calibration with two sources of natural gas, and it assumes that any variation in the composition of natural gas is due to a combination of these two sources. It is therefore convenient to select an “H” gas and an “L” gas for this purpose. If only these two sources (and hydrogen) are randomly blended, this method would exhibit virtually no detection error. The final estimator whose performance will be evaluated is based on the combined measurement of the three mentioned variables (thermal conductivity, speed of sound and density). This estimator was calibrated with the three types of natural gas whose thermal conductivity is shown in Figure 2 (an “H” gas, an “L” gas, and a type of biomethane). Its underlying assumption is that any variation in the constitution of natural gas is due to a ternary mixture of these three sorts. This estimator is therefore expected to be much more robust against composition variations, detecting both the concentration of hydrogen and the combustion enthalpy with a high precision—in spite of there being seven different sources of natural gas, mixed binarily in random proportions. The accuracy of these four estimation methods was evaluated with the root-mean-square error (Equation 15) in the prediction of the variables of interest (combustion enthalpy and molar fraction of hydrogen) for all the randomly generated mixtures. ϵ R M S = 1 N ∑ x p r e d i c t − x t r u e x t r u e 2 (15)

The previous comparison assumes that the measurement of the physical variables is free from experimental error, which is unrealistic. For this reason, a second study was performed, in which the outcome of each physical measurement is distorted by a Gaussian distribution with varying degrees of standard deviation (0.5%, 1.0% and 2.0% of the mean). Therefore, for every randomly generated gas mixture, a series of 30 “measurements” of the physical variables under the influence of said cases of experimental error were performed, in order to assess the robustness of the proposed detection methods in a realistic scenario. It is expected that the three-variable estimator will suffer a higher loss of accuracy under the effect of experimental uncertainty than the two- or one-variable ones.

We proposed in Section 2.2 to simplify the calculation of the thermal conductivity of a mixture (Equation 7) as a quadratic polynomial (Equation 9). It can be seen in Figure 2 that the coefficients of such a polynomial can be fitted with a least-squares procedure to predict the thermal conductivity with high precision. These coefficients were fitted using the simulated curves of thermal conductivity vs amount fraction of hydrogen for three selected sources of natural gas (one ‘H’-type, one ‘L’-type and a type of biomethane). Such a polynomial fit could be obtained with experimental data by measuring the thermal conductivity of each independent natural gas source upon hydrogen addition in the concentration range of interest. This quadratic model facilitates the solution of a system of equations to a great extent. Some have proposed to use a linear fit for the thermal conductivity (Huber et al., 2018), but this would result in a loss of precision.

In order to find an appropriate fit of the constants in Equation 9, we suggest using an optimisation routine with a first guess based on the separate regressions of the conductivity-concentration curves of the three independent natural gas sources used for calibration. The substitution for each separate combination of natural gas and hydrogen leads in each case to a simplification of the complete polynomial, as shown in Equation 16. The first case corresponds to x 2 = 0 ; x 1 = 1 − x H ; the second case corresponds to x 1 = 0 ; x 2 = 1 − x H ; and the third case to x 1 = 0 ; x 2 = 0 . From each case, a set of three constants can be directly obtained by a standard polynomial regression; however, the nine constants that are obtained cannot be used directly to calculate the ten constants required by Equation 9. Our suggestion is to set all cross-coefficients κ 110 , κ 101 , κ 011 to zero, which enables obtaining the initial guess of the rest of the constants by direct substitution. The actual coefficients can be then obtained by a minimisation of the sum of the root-mean square error of all three curves against the complete polynomial. κ 200 + κ 020 − κ 110 ⁢ x H 2 + κ 100 − 2 κ 020 + κ 110 − κ 010 ⁢ x H + κ 020 + κ 010 + κ 000 = κ κ 200 + κ 002 − κ 101 ⁢ x H 2 + κ 100 − 2 κ 002 + κ 101 − κ 001 ⁢ x H + κ 002 + κ 001 + κ 000 = κ κ 200 ⁢ x H 2 + κ 100 ⁢ x H + κ 000 = κ (16)

The variation in the speed of sound, density, thermal conductivity and viscosity of all the randomly generated natural gas compositions can be observed in Figure 3. One can appreciate how the variation in the composition of natural gas has a non-negligible impact in the speed of sound and density, a mild effect in the thermal conductivity, and a dominant role in the viscosity. It is noticeable that, for instance, a reading of 460 m/s in the speed of sound corresponds to 10 mol% H 2 in one end and to 20 mol% H 2 in the opposite end of the spectrum of natural gas composition. If such a variation in the constitution of natural gas is present in a real application, the measurement of just the speed of sound or density would be inadequate to predict the amount fraction of hydrogen. On the other hand, the high sensitivity of the thermal conductivity to the presence of hydrogen, combined with a relative insensitivity to the variation in the composition of natural gas, makes this variable very appealing for estimating the hydrogen content. The opposite case is observed in the viscosity, whose variation with the amount fraction of hydrogen (at least in the range until 30 mol%) is rather negligible, whereas its sensitivity to the composition of natural gas is very high.

Having calculated the thermodynamic properties of several random gas mixtures, the accuracy of the methods proposed in Section 2.3. Can be evaluated. In spite of ignoring the variations in the composition of natural gas, the one-variable estimator based on the thermal conductivity predicts the hydrogen content relatively well, with an expected deviation around ±2 mol%. Such a deviation would be inadequate for low H 2 concentrations (e.g., 5 mol%) but uncritical for higher concentrations (e.g., 30 mol%). Figure 4 shows the root mean square error of the estimator as a function of the hydrogen concentration. The expected error in the estimation of the combustion enthalpy is close to 5%, which could be tolerable for household billing applications. By measuring an additional variable to the thermal conductivity (whether density or the speed of sound), the two-variable estimator was able to reduce the error in the prediction of the amount fraction of hydrogen significantly. Nevertheless, there seems to be no improvement in the estimation of the combustion enthalpy. The three-variable estimator, whilst retaining a comparable error in the prediction of the amount fraction of hydrogen in comparison to the two-variable cases, was able to further improve the accuracy of the prediction of the combustion enthalpy.

If the effect of a random measurement uncertainty is considered, the simultaneous measurement of more variables does not automatically lead to a more robust estimation. As Figure 5 shows, the impact of the simulated noise levels is negligible for the estimator based on only one variable (thermal conductivity). A different situation is observed for the estimators based on two variables, as shown in Figures 6, 7. The estimation of the hydrogen content was increasingly distorted with the introduced noise, whereas the estimation of the combustion enthalpy remained unaffected. With a measurement error of 2% in both variables, the estimator based on speed of sound and thermal conductivity performed no better than the one based on thermal conductivity alone, and the one based on thermal conductivity and density performed only slightly better. The final case of a three-variable estimator showed, as expected, a very high sensitivity to the introduced noise (Figures 8). Only if the noise in all three variables is kept below 1% does it show an advantage compared to the one-variable estimator.

Considering all these effects, the combination of thermal conductivity and either density or speed of sound seem to offer a good trade-off between robustness against random variations and robustness against noise. These two estimators outperform the measurement of thermal conductivity alone, as long as the total measurement uncertainty of both variables is kept around or below 1%. Although the measurement of thermal conductivity and density seems to be less sensitive to measurement noise with respect to the case with the speed of sound, the measurement of density is affected by the propagated error proceeding from the temperature and pressure compensation that it requires (whereas the speed of sound only requires a temperature compensation). A further advantage of measuring the speed of sound is that, if it is determined by an ultrasonic system, such a device will naturally also provide the measurement of the flow velocity, which is necessary for computing the energy consumption (Equation 3).

It has been suggested that the heating value of natural gas (without blending hydrogen into it) can be correlated to some thermodynamic properties (Huber et al., 2018). In particular, the higher heating value ( H H V = − ρ M Δ ° H c ) is said to correlate with density, whereas the Wobbe index ( W o = H H V / ρ / ρ a i r ) with viscosity. With the randomly calculated mixtures in our Monte-Carlo simulation, we do not observe a clear tendency between these variables, as Figure 9 shows. However, we cannot discard the possibility of a correlation based on a more sophisticated combination of these variables and perhaps some additional ones.