Response time correction of slow response sensor data by deconvolution of the growth-law equation

,


Introduction
High resolution in situ data are crucial to observe high variability in environmental processes when surrounding environmental parameters are continuously changing.Many contemporary measurement techniques have a limited response time due to signal convolution inherited from a diffusion process, such as in Vaisala radiosondes (Miloshevich et al., 2004) or in continuous flow analysis of ice cores (Faïn et al., 2014).In oceanic sciences, measurement of dissolved analytes often requires an extraction technique based on membrane separation, where the property of interest (a solute) equilibrate across a membrane separating the medium of interest (a solvent) from the medium where the actual measurement takes place (a measurement chamber).This makes the sensor response time (RT) directly governed by how fast the analyte of interest can diffuse through the membrane.
This process is mainly driven by the difference in partial pressure between the two media and can be relatively slow.This results in high sensor RTs, leading to unwanted spatial and temporal ambiguities in recorded signals for sensors used in profiling (Miloshevich et al., 2004), on moving platforms (Bittig et al., 2014;Canning et al., 2021) or deployed in dynamic environments (Atamanchuk et al., 2015).Herein, we refer to sensors with this particular design as Equilibrium Based (EB) sensors and we seek to establish a robust, simple and predictable method for correcting high RT induced errors in data from these sensors.
Considering an EB sensor during operation, we define u a (t) as the instantaneous ambient partial pressure of interest and u m (t) as the partial pressure within the measuring chamber of an EB sensor, where the measurement occurs.In this situation, a model of u m (t) as a function of time can be obtained via the growth-law equation (Miloshevich et al., 2004), which describes diffusion of the property of u a (t) through the separating barrier (in this case the membrane): where k is a sensor specific growth coefficient, which determines how fast a change in u a (t) will be reflected in u m (t).The RT of EB sensors are often given in τ 63 = 1/k, which corresponds to the time the sensor requires to achieve 63% (one e-folding) of an instantaneous step-change in ambient concentration.If k in Eq. 1 is sufficiently small (i.e.τ 63 is large), the diffusion will be slow and any fast fluctuations in u a (t), will be smeared out in time.
A numerical technique has already been proposed to recover fast fluctuations in u a (t) from measurements of u m (t) using a closed form piece-wise solution to Eq. 1 (Miloshevich et al., 2004).However, due to the ill-posed nature (see e.g.Tikhonov et al., 1977) of the forward model, errors in the measurements will be amplified when reconstructing u a (t).Miloshevich et al., (2004) counteracts this using an iterative algorithm that minimizes third derivatives to obtain locally smooth (noise-free) time-series prior to the reconstruction of u a (t).While this and similar methods seems to usually work well in practice (Bittig et al., 2014;Canning et al., 2021;Fiedler et al., 2013;Miloshevich et al., 2004), it is difficult to determine the uncertainty of the estimate, as the iterative scheme does not model error behavior.Predicting the expected solution of the iterative estimator is also difficult and there is no straightforward way of choosing suitable smoothing parameters.These are important attributes for the reliability of solutions to these types of problems, due to the error amplification that occurs during deconvolution.
Herein, we establish an alternative method for estimating u a (t) from a measurement of u m (t).This solution is based on the framework of statistical inverse problems and linear regression.Using a weighted linear least squares estimator, the growth law equation as measurement model, and a sparsity regularized solution, we are able to take into account uncertainties in the measurements, provide an intuitive and/or automated way of specifying an a priori assumption for the expected solution, and determine the uncertainty of the estimate.This approach also enables us to evaluate the self-consistency of the solution and detect potential instrument ::::: model : and/or measurement issues.A time-dependent k can also be employed, which suits membranes with varying permeability (e.g.where k is a function of temperature).We show that automated L-curve analysis produces well regularized solutions, thereby reducing the number of input parameters to sensor response time, measurement uncertainty, and the measurements themselves.The robustness/functionality of our technique was validated using simulated data, laboratory experiment data, and comparison of simultaneous field data from a prototype fast response Diffusion Rate Based (DRB) sensor (Grilli et al., 2018) and a conventional slow response EB sensor in a challenging Arctic environment.

Method
We assume that the relationship between observed quantity u a (t) and measured quantity u m (t) are governed by the growthlaw equation as given in Eq. 1. Estimating u a (t) from u m (t) is an inverse problem (Kaipio and Somersalo, 2006;Aster et al., 2019)meaning that , ::: and : a small uncertainty in u m (t) will result in a much larger uncertainty in the estimate of u a (t), making it impossible to obtain accurate estimates of u a (t) without prior assumptions.
To formulate the measurement equation (Eq. 1) as an inverse problem that can be solved numerically, we need to discretize the theory, model the uncertainty of the measurements, and establish a means for regularizing the solution by assuming some level of smoothness :::: (prior ::::::::::: assumptions).We will denote estimates of u a (t) and u m (t) as ûa (t) and ûm (t).Measurements of u m (t) will be noted as m(t).Each element of the following steps is illustrated in Figure 1.
We discretize Eq. 1, using a time-symmetric numerical derivative operator ::: (e.g.::::::: Chung, :::: 2010 : ): We have used the following short-hand to simplify notation: u m (t i ) = u i and u a (t i ) = a i .Here t i is an evenly sampled grid of times and ∆t i = t i+1 − t i is the sample spacing.We refer to t i as model time : , and for simplicity , :: we : assume that this is on a regular grid ∆t i = ∆t with a constant time step.Note that the growth coefficient k i = k(t i ) can vary as a function of time.
We assume that sensor measurements m j of the quantity u m (t) obtained at times t j (see Figure 1) have additive independently distributed zero-mean Gaussian random noise: where ξ j ∼ N (0, σ 2 j ) with σ 2 j the variance of each measurement, which in practical applications can be estimated directly from the data or by using the : a known sensor accuracy.The t j is the measurement time, which refers to the points in time where measurements are obtained.We obtain ::::::: Coupling ::::::: between : u m (t i ) ::: and :::::::::::: measurements :::: m j :: is ::::::: obtained : through gridded re-samplingof m j (see Figure 1).Note that the measurement time-steps t j do not need to be regularly spaced, nor coincide with the model times t i .
To reliably estimate u a (t), we need to regularize the solution by assuming some kind of smoothness for this function.A common a priori assumption in this situation is to assume small second derivatives of u a (t), corresponding to the second order Tikhonov regularization scheme (Tikhonov and Arsenin, 1977).Although this provides acceptable solutions, the choice of the regularization parameter (i.e.adjusting the amount of regularization applied) is not particularly intuitive.Since our method is intended for a variety of domains where validation can be challenging, we have chosen to employ a different regularization method, where the regularization parameter relates directly to simple, real world characteristics and the ability of the instrument to resolve the ambient environment.
Model sparsity regularization (see e.g.Hastie et al., 2015) provides an intuitive model regularization by assuming that the observed quantity can be thoroughly explained by a reduced number of samples in some domain.In our case, we have used time domain sparsity, which translates to setting the number of model time (t i ) steps N smaller than the number of measurement time (t j ) steps M. The a priori assumption we make to achieve this is that the observed quantity can only change with a time and change piece-wise linearly between these points.Using this approach, an optimal regularization parameter becomes the lowest ∆t at which the observed quantity can change significantly.This means that choosing the regularization parameter can be done based on domain specific knowledge or scientific requirements for temporal resolution, within the limitations posed by the ill-posed nature of the problem and sensor performance.
We can now express the theory, relationship between the measurements and the theory, and the smoothness assumption in matrix form as follows: Eq. 2).N is the number of model grid points i and M is the number of measurements m j .A ∈ R N ×N and U ∈ R N ×N express the growth-law relationship between discretized a i and u i as given in Eq. 2. The matrix V ∈ R M ×N handles the regularization and the relationship between measurements and the discretized model of u m (t).The constant γ σ −1 j is a numerically large weighting constant ::: (we :::: used ::::::::::::::: γ = 2∆t • 10 5 σ −1 ), which ensures that when solving this linear equation, the solution of the growth-law equation will have more weight than the measurements.In other words, the solution satisfies the growth-law equation nearly exactly, while the measurements are allowed to deviate from the model according to measurement uncertainty.The definitions of the matrices are as follows, where k i is the growth coefficient: In matrix U we also need to express the time derivative of u m (t) and consider edge effects: The matrix V ∈ R M ×N relates concentration u i to measurements of this concentration m j (see Figure 1 and Eq. 2 and 3).We use a weighted linear interpolation between grid points t i when assigning measurements to the model: Where ∆t is the model timestep and regularization parameter (see Eq. 4).It is now possible to obtain a maximum a posteriori estimate of u a (t) and u m (t) by solving for the least-squares solution to matrix Eq. 5: The vector x = [â; û] contains the maximum a posteriori estimate of vectors a and u.The matrix G is described in Eq.
5-9.The estimate â of u a (t) is the primary interest; however, the solution also produces an estimate û of u m (t).This can be a useful side-product for detecting outliers from fit residuals or other issues with the measurements (as shown in the field experiment).
Equation 10 allows the estimate to also be negative, which can be unwanted if the observed quantity is known positive.In such cases, it is possible to apply a non-negativity constraint using a non-negative least-squares solver (Lawson and Hanson, 1995).
As the uncertainties of measurements are already included in the theory matrix G, the a posteriori uncertainty of the solution is contained in the covariance matrix Σ MAP , which can be obtained as follows: This uncertainty includes the prior assumption of smoothness.

Simulation and ∆t determination
To test the numerical validity of our method and develop a regularization parameter selection tool, we used a toy model.This gives us the possibility to prove that the method gives well behaved consistent solutions as we know the correct results and control all input variables.
We defined the simulated concentration u a (t) (see also Figure 1) as a step-wise change in partial pressure: where units for time and partial pressure are arbitrary.While this is not a realistic scenario encountered under field conditions, step-change simulations is a conventional calibration.It is also the most challenging scenario for testing our method, since it directly violates our smoothness assumption.
As A, U , and m of the matrix equation (Eq.6) are now known, only the regularization parameter ∆t in V (Eq.9) needs to be defined to obtain the gridded estimate of time and RT-corrected measurements(i.e.û and â) for the simulated model.
For solutions regularized through smoothing, a well regularized solution captures a balance between smoothness and model fit residuals.Or in more practical terms, the provided solutions are sharp enough so critical information is not lost (i.e., detection of short term signal fluctuations are not removed by smoothing integration), but with reasonable noise and uncertainty estimates.
We have chosen a heuristic approach to this optimization problem, by applying the L-curve criterion (see Hansen, 2001) Statistical methods based on Bayesian probability (Ando, 2010), such as the Bayesian information criterion were also tested with similar results.We chose to apply the L-curve criterion due to its robustness and ability to intuitively display the effect the regularization parameter has on the solution, which we believe is an advantage in practical applications of our method.
In our case, the L-curve criterion involves calculating a norm E s , which measures how noisy the estimate is and a norm E m , which measures how large the fit residual is (i.e., an estimate of how well the model describes the measurements).These norms are calculated for a set of different regularization parameters, which are compared in a log-log plot (see Figure 2a) where the data points align to trace a curve that resembles the letter "L".The under-regularized (or too noisy) solutions are found in the upper left corner where perturbation errors dominate.The over-regularized (or over-smoothed) solutions are located in the lower right corner, where regularization errors dominate.Good regularization parameters are located in the middle of the bend or kink of the L, where smoothness and sharpness are well balanced, limiting both noise and fit residuals.
We have used the first-order differences of the maximum a posteriori solution â as the norm measuring solution noise: and to approximate the fit residual norm, we use: where m j is the measurement of the quantity u m (t): which comes directly from the least squares solution (Eq.10) of matrix equation 6, where corresponds to V ji (see Eq. 9) but without scaling for measurement error standard deviation.In essence, ûm (t j ) is the best fit model for the measured quantity u m (t) is obtained using linear interpolation in time from the least squares estimates in vector the ::::: bend :: of ::: the ::: L).
Using ∆t =1.35, we can also inspect how a well-regularized solution is affected by edges and missing measurements (shown in Figure 3d).The missing measurements resulted in an increased uncertainty in the region where there are no measurements, but the maximum a posteriori estimate still provided a reasonable solution.Uncertainties grow near the edges of the measurement as expected, since there are no measurements before or after the edges which contain information about u a (t).

Laboratory experiment
We evaluated our proposed technique in a controlled laboratory experiment using a Contros HydroC CH 4 EB methane sensor by exposing the sensor to step changes (similar to the numerical experiment) in concentration.We connected the instrument to an air tight water tank (12.9 L) with a small headspace (∼0.25L) via hoses where water was pumped at 6.25 L min −1 and kept at constant temperature (22 o C) (see setup in Appendix B).The setup first ran for two hours to ensure stable temperature and flow.In this period, the sampling rate was 60 seconds, while it was changed to 2 seconds for the rest of the experiment to provide high measurement resolution for the step changes.
Four step-changes (two up and two down) were approximated by opening the lid and adding either 0.2 L of methane enriched or ultrapure water.The RT of the sensor was determined to 23 minutes at 22 o C (k=0.000725s −1 ) prior to the experiment following standard calibration procedure of the sensor, and parameters affecting membrane permeation was controlled (i.e. water flow rate over membrane surface and water temperature).We therefore used the signal noise calculated ::::: signal ::::: noise ::::::::: represented : by the standard deviation of the single point finite differences as measurement uncertainty . :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: Input uncertainty was lower for the first 2 hours as the lower sampling rate reduces signal noise due to a longer internal averaging period.∆t was determined to 179 s using the automatic ∆t selection based on the L-curve criterion (see Appendix  A).This also gave well confined fit residuals, with a difference in signal noise at around 7000 s due to the change in input noise (caused by changes in sampling rate from 60 seconds to 2 seconds after 2 hours, Figure 4a and b).
At the first step change, the RT-corrected concentration rapidly increased from around 2.6 µAtm to ∼41 µAtm (Figure 4c and d).This was followed by a slower increase taking place over around 30 minutes up to ∼47 µAtm and then a slow decrease for another 30 minutes down to ∼45 before the next step change.The following step increase and subsequent two step decreases followed the same pattern, which can be explained by three processes indicated in Figure 4d): The first rapid increase (process 1) results from the initial turbulent mixing caused by the abrupt addition of methane enriched water to the tank.This is followed by a slower diffusive mixing phase occurring after the water has settled (process 2).While these processes are occurring, there is also a gradual diffusion of methane to the headspace, which is shown as decreasing concentrations in approximately the last half of the plateau periods (process 3).As expected, this decrease was faster for higher concentrations, due to the larger concentration gradients across the water/headspace interface.The step decreases show the same behavior, although with process 2 inverted.
Estimated uncertainty of the RT-corrected data averaged 0.64 µAtm (95% confidence) which is roughly double the raw data noise of 0.29 µAtm.Due to the long concentration plateaus, the balanced ∆t lies in a quite strongly regularized solution.
Nonetheless, the de-convolved instrument data gives a considerably better representation of the step-changes with a relatively small uncertainty estimate and reveals known features of the experiment setup (processes 1-3 in Figure 4d) which is obscured in the convolved data.

Field experiment
Continuted evaluation of our proposed technique was applied under more challenging conditions in a field based study using simultaneous data from two different methane sensors towed over an intense seabed methane seep site offshore West Spitsbergen (Jansson et al., 2019).A slow response EB Contros HISEM CH 4 and a fast response Membrane Inlet Laser Spectrometer (MILS) DRB sensor (Grilli et al., 2018) were mounted on a metal frame and dragged at various heights over the seabed (∼20-300 m) at 0.4-1.1 m s −1 in an area with many hydro-acoustically mapped methane seeps.The rapid and large variability in methane concentration, direct comparison with the DRB sensor, and the particularly high RT of the EB sensor in cold water made this an ideal test scenario for field based applications.
The measurement uncertainty (ξ j :: σ j in Eq. 14) was set to either the estimated raw data noise or to the stated sensor accuracy after equilibrium is achieved, depending on which of these parameters was higher.The EB accuracy is stated to 3% using the ISO 5725-1 definition of accuracy, which involves both random and systematic errors.High concentrations during the field experiment made the 3% sensor accuracy our main input parameter for measurement uncertainty.

∆t determination
We produced an L-curve from a set of estimates of u a (t) with ∆ts ranging from 10 to 550 seconds (Figure 5a).Using a polynomial spline to approximate the maximum curvature point (see Appendix A) we found the optimal ∆t to be 55 s, which is in good agreement with our visual inspection of the L-curve.Upon inspection of the fit residual plot (Figure 5b), we observed large spikes at several time points, meaning that the model failed to describe the transformation between the measurements m j and the RT-corrected estimate u a (t i ).Inspection of raw data (red line in Figure 5b) uncovered sharp signatures in the measured dissolved concentration at these instances -too sharp to be a real signal (as these should have been convolved by the instrumental function).We concluded that there was an unidentified problem with the EB sensor system at these instances, most likely related to power draw, pump failure or other instrumental artifacts.Removing these problematic sections and redoing the estimates provided an L-curve with an unchanged maximum curvature location.We therefore kept using ∆t=55 s in the final deconvolution which now gave a solution with approximately Gaussian distributed fit residuals without spikes (Figure 5d), meaning that we have a self-consistent and valid solution.
The high R 2 of the RT-corrected data confirms that the EB sensor captured most of the variability in dissolved methane; however, there is a slight bias in the differences between the two data sets (α=0.82).Inspecting the absolute differences reveals that the RT-corrected EB data have more moderate concentrations during periods of very strong variability.This can partly be explained by the inherent smoothing of a sparse model (for larger ∆t).In theory, increasing the model complexity  (reducing ∆t) should return a slope closer to 1 and reduce differences.However, our attempt at decreasing ∆t for achieving this did not improve the slope, but increased the noise as expected.The flat slope could also at least partly originate from the previously problematic sections (spikes in Figure 5) in the EB sensor data (arrows in Figure 5b).Even though we ignored the data at these intervals, the offset in absolute concentration still affects our end-estimate.Another explanation could be an overestimation of the k of the EB sensor in the laboratory procedure prior to the field campaign.Indeed, when using a lower k, e.g.corresponding to τ 63 =3600 s (keeping ∆t=55 s), which matches better with calibration results for similar sensors, the slope goes to 1 and differences are evenly distributed between high and low methane concentrations.The small slope offset could also be a combined effect of the above iterated reasons.
Using the framework of inverse theory allows us to model error behavior, enabling a comparison of the uncertainty estimates of the two sensors (shaded regions in Figure 6).For the DRB sensor data, we used the stated 12% accuracy as the uncertainty 310 estimate (Grilli et al., 2018).For the RT-corrected EB data, we used the 95% confidence of the deconvoluted estimate using input measurement uncertainty, resulting in a median uncertainty range of 22%.Linearly interpolating the RT-corrected EB data onto DRB data time and taking the mutual error bounds into account, the two data sets agree within the uncertainties 92% of the time.This is despite the lower resolution of the RT-corrected data and other error sources (some of them described above) and we consider this a successful result.
The DRB data has a lower median relative (%) uncertainty estimate, but to compare these relative uncertainty estimates directly can be slightly misleading as the relative uncertainty estimate of the RT-corrected data varies in time (Figure 6d).This is due to the EB sensor convolution occurs prior to the time when the actual measurement (including measurement uncertainty) takes place (see Figure 1b).Since input measurement uncertainty mostly follows 3% of measured (already convolved) value, the uncertainty estimate becomes a function of the EB raw data.Consequently, due to the raw data hysteresis, the uncertainty becomes lower for increasing concentrations compared to decreasing concentrations, and vice versa.We can observe this at ∼00:15 and ∼07:15, when the RT-corrected concentration increases dramatically, but the uncertainty estimate is still relatively small (Figure 6a and d).On the other hand, between 04:40 and 05:30, the RT-corrected concentration data is relatively low and constant, but the uncertainty estimate is large and shrinks slowly due to the slow decrease in u m (t).Comparing the error bounds of the data-sets using the relative uncertainty is therefore a simplification because of the raw-data -inherited RTcorrected uncertainty estimate.Overall, this adversely affects the relative uncertainty estimate of the RT-corrected data set, since high error bounds inherited from the raw-data hysteresis during decreasing concentrations is divided by RT-corrected low concentration values which results in high relative uncertainties.A more narrowly defined (if possible) input uncertainty estimate for the EB sensor could help constrain this uncertainty.

Conclusions
We presented and successfully applied a new RT-correction algorithm for membrane based sensors through a deconvolution of the growth-law equation using the framework of statistical inverse problems.The method requires few and well-defined input parameters, allows the user to identify measurement issues, models error propagation and uses a regularization parameter which relates directly to the resolution of the response time corrected data.Functionality testing was done using both a laboratory and a field experiment.Results from the laboratory experiment uncovered features of the experimental setup which were obscured by convolution in the raw data and the field experiment demonstrated the robustness of the algorithm under challenging environmental conditions.In both tests, the sensors ability to describe rapid variability was significantly improved and better constraints on input uncertainty and response time are areas which can potentially further enhance results.
This method and validation experiments using the Contros/HISEM sensors uncovers a new set of applications for these and similar sensors, such as ship-based profiling/towing and monitoring highly dynamic domains.Conventional EB sensors are also more abundant and affordable compared to more specialized equipment, increasing the availability and possibilities for scientists requiring high resolution data to solve their research questions.Additionally, we believe this deconvolution method could be applicable to other measurement techniques as well, where diffusion processes hampers response time.

Figure 2 .
Figure 2. a) L-curve for sweep of different ∆t values for estimated property (RT correction) of simulated measurements given by mj = um(tj) + ξj.The y-axis, Es is the noise in the data given by the step difference between adjacent data-points which is high for models that are too complex (lower ∆t/higher number of data points) in the model.The x-axis is the fit error residual Em, which shows how well the results (the uas) explains the measurements m.The latter is high for too sparse models (high ∆t/low number of data points).b), c), and d) show fit residuals û(t j ) − mj for each point in measurement time using ∆t=0.25 (b), ∆t=5 (c), and ∆t=1.35.

Figure 3 .
Figure 3.Estimated property using simulated data and different regularization parameters (time intervals).The simulated measurements mj = um(tj) + ξj are shown with red points, and the estimated property ûa(t) is shown with a blue line with 2-σ uncertainty indicated with a light blue line.Panel a) shows a high temporal resolution estimate, which is very noisy, due to the ill-posed nature of the deconvolution problem, b) shows a too low temporal resolution estimate, which fails to capture the sharp transition at t = 5, and c) shows a good balance between estimation errors and time resolution.Panel d) shows estimated property for a simulated data set with missing measurements between t = 15-17, which results in an increased error around the missing measurements for the estimated quantity.Error estimates are given as 95% confidence intervals.

Figure 4 .
Figure 4. a) L-curve for sweep of different ∆t values for estimated RT corrected EB sensor data.The y-axis, Es, is the noise in the datagiven by the step difference between adjacent data-points which is high for models that are too complex (lower ∆t/higher number of data points) in the model.The x-axis is the fit error residual, ûm(t j ) − mj, which shows how well the best fit model explains the measurements.b) Fit residuals ûm(t j ) − mj (black line) for each point in measurement time using ∆t=179 s. c) show the result of the deconvolution (black line), uncertainty estimate (grey shading) and raw EB sensor data (blue line).d) is a zoomed in version of c) with area specified in c).

Figure 5 .
Figure 5. a) L-curve for sweep of different ∆t values for estimated RT corrected EB sensor data.The y-axis, Es, is the noise in the datagiven by the step difference between adjacent data-points which is high for models that are too complex (lower ∆t/higher number of data points) in the model.The x-axis is the fit error residual, ûm(t j ) − mj, which shows how well the best fit model explains the measurements.b) Fit residuals ûm(t j ) − mj (blue line) for each point in measurement time using ∆t=55 s and raw EB sensor data (red line).Residual spikes attributed to a malfunctioning of the EB sensor are indicated by the black arrows (and zoom in oval inlet) c) and d) reported the results from the same analysis as a and b, but on the dataset where the problematic regions defined in figure b were removed.

Figure 6 .
Figure 6.a) Field data from the DRB (yellow), EB (red), and RT-Corrected EB (blue) sensors.b) and c) show direct comparison between DRB and the EB sensor data.Coefficient of determination (R 2 ), Mean Absolute Error (MAE), and slope angle (α) is given for comparison between the DRB sensor data and either raw (b) or RT-corrected (c) EB sensor data.d) show the error estimate of the RT-corrected signal, ua(ti).

Figure
Figure B1.a) Schematic representation of the experiment setup and b) picture of the experiment setup.The tank had a air tight dome shaped lid with a small headspace (∼0.3 L) and was 27.5 cm high and 24.5 cm diameter.Room temperature was controlled and kept constant.