Extreme Indian rainfall under El Niño]More extreme Indian monsoon daily rainfall in El Niño summers

[1]Spencer A<EMAIL_ADDRESS>2]Destiny Zamir<EMAIL_ADDRESS>3,4]Adam H<EMAIL_ADDRESS>3]Michela<EMAIL_ADDRESS>3]Mark A<EMAIL_ADDRESS>4]Michael K<EMAIL_ADDRESS>5]Fiaz<EMAIL_ADDRESS>
*[1]Department of Earth and Atmospheric Sciences, City College of New York, 160 Convent Avenue, New York, 10031, New York, United States

[2]Department of Civil Engineering & Engineering Mechanics, Columbia University, 500 W. 120th Street #610, New York, 10027, New York, United States

[3]Lamont-Doherty Earth Observatory, Columbia University, 61 Route 9W, Palisades, 10964, New York, United States

[4]Department of Applied Physics and Applied Mathematics, Columbia University, 500 W. 120th Street #200, New York, 10027, New York, United States

[5]Department of Atmospheric and Oceanic Sciences, UCLA,  Los Angeles, 90095, California, United States

Extreme rainfall in the Indian summer monsoon can be destructive and deadly <cit.>.  Although / events in the equatorial Pacific make dry days and whole summers more likely throughout India <cit.>, their influence on daily extremes is not well established.  Despite this summer-mean drying effect, we show using observational data spanning 1901-2020 that / increases extreme rainfall likelihoods within monsoonal India, especially in the the summer's core rainy areas of central-eastern India and the narrow southwestern coastal band.  Conversely, extremes are broadly suppressed in the drier southeast and far northwest, and more moderate accumulations are inhibited throughout the domain.  These rainfall signals appear driven by corresponding ones in convective buoyancy, provided both the undilute instability of near-surface air and its dilution by mixing with drier air above are accounted for <cit.>.  When the summer ENSO state is predicted from a seasonal forecast ensemble initialized in May, the extreme rainfall patterns broadly persist, suggesting the potential for skillful seasonal forecasts.  The framework of analyzing the full distributions of rainfall and convective buoyancy could be usefully applied to hourly extremes, other tropical regions under ENSO, other variability modes, and to trends in extreme rainfall under climate change.

The influence of the / (ENSO) on total summer rainfall in India has been recognized since the early 20th century <cit.>.  Anomalously warm equatorial Pacific sea surface temperatures (SSTs) in / summers generate anomalous season-mean ascent locally but compensating descent over much of the rest of the tropics including India.

One might expect summer-mean subsidence to inhibit moist convection of all intensities, extremes included, reducing both the mean and variance of daily rainfall.  Indeed, recent studies <cit.> establish that / reduces the average accumulation on rainy days—the rainfall intensity—within southeastern India and in Rajasthan state in the northwest.  These are the summer monsoon's climatologically dry areas, and within them during / summers it tends to rain even less overall, less often, and less intensely  <cit.>.

But the opposite rainfall intensity signal emerges in the summer monsoon's climatologically rainy areas <cit.>: the broad Central Monsoon Zone in central India and the narrow southwestern coastal band of peninsular India.  Within them during / summers it tends to rain less overall and less often—but more intensely when the rain comes.

Whether these signals in rainfall intensity stem from destructive extremes <cit.> or beneficial moderate events <cit.>—and mechanistically what causes them—remain unanswered.  Observational studies have investigated ENSO's influence on daily or hourly rainfall extremes in other tropical domains <cit.> but not the Indian summer monsoon.

§ ENSO INFLUENCES ON EXTREME DAILY RAINFALL

Our primary metric of extreme daily rainfall is the cutoff accumulation <cit.>, computed as the ratio of the variance to the mean of daily rainfall across summer days (Methods).  This metric is motivated by the fact that daily rainfall reliably follows a gamma distribution <cit.>: probabilities follow a shallow power-law below the cutoff before dropping off exponentially above it (Fig. <ref>a).  Increasing the cutoff, which in the exact gamma limit is identical to the scale parameter, therefore makes all extreme rain rates more likely.  Our primary dataset is the Indian Meteorological Department daily 0.25\(\times\)0.25 gridded rainfall product spanning 1901-2020 derived from a dense, temporally varying in situ rain gauge network <cit.>.  We use all 122 days within each June-July-August-September (JJAS) summer season and all gridpoints within the “monsoonal India” domain <cit.> (Methods).  As shown below, key results are robust across alternative metrics of extremes and to another daily rainfall dataset.

The climatological cutoff spans from 11.8 to 75.9 mm across all monsoonal India gridpoints (Fig. <ref>b), with a spatial structure that largely tracks that of summer-mean rainfall: high values within the southwestern coastal band, a sharp gradient moving east across the Western Ghats mountain range to much smaller values in the southeast, intermediate values north thereof in most of the Central Monsoon Zone, and low values again in the far northwest.  This pattern emerges in many properties of summer monsoon rainfall from daily to interannual timescales <cit.> (Fig. <ref>a; overlaid contour in Fig. <ref>b; Extended Data Fig. <ref>).  The cutoff and other metrics of extremes are large compared to the mean in Gujarat (along the Arabian Sea coast within the Central Monsoon Zone), which is rather dry most summers but in others receives synoptically driven extreme accumulations <cit.> (Extended Data Fig. <ref>).

(a) Daily rainfall empirical probability distributions (Methods), in log-log space, combining all June-July-August-September (JJAS) days 1901-2020 and all gridpoints within all of monsoonal India (black), the Central Monsoon Zone (blue), the Western Ghats coastal band (green), or Southeast India (brown).  Overlaid circles indicate the climatological cutoff.
(b) Local climatological cutoff daily rainfall (mm) across all JJAS days 1901-2020 according to the colorbar.
The portions of northern and northeast India shaded grey are excluded from all analyses (Methods).
The dotted black contour denotes the median across gridpoints of the climatological 99th percentile daily rain rate (63.1 /).
(c) Pearson correlation coefficient between the JJAS NINO3.4 value and each gridpoint's cutoff rain rate.  The dotted black contour denotes the JJAS climatological rainfall median (5.9 /) across gridpoints.  Semi-transparent shading indicates statistical insignificance at the 95% confidence level based on a block-bootstrap method.  All timeseries are linearly detrended over 1901-2020 prior to correlations being computed.

Based on lag-zero correlations of each gridpoint's JJAS cutoff with the standard NINO3.4 index of ENSO (Methods), extreme daily rainfall becomes more likely over much of the southwestern coastal band and central India as NINO3.4 increases but less likely over much of the southeast and far northwest (Fig. <ref>c).  This pattern resembles the climatological summer-mean rainfall (overlaid contour in Fig. <ref>c; Extended Data Fig. <ref>) and the established <cit.> ENSO-driven patterns in rainfall intensity.  As one quantification of this resemblance, of the 581 gridpoints with statistically significant positive correlations, 75% occur in gridpoints with climatological summer-mean rainfall above the median across gridpoints, and of the 384 gridpoints with significant negative correlations, 99% occur in gridpoints below it.  As another, linear regressions of each gridpoint's NINO3.4-cutoff correlation / its summer-mean climatological rainfall are significantly positive for the whole domain and within each subregion (Extended Data Fig. <ref>).

For each of the four subregions of India outlined in Fig. <ref>, we compute an aggregate cutoff by combining all gridpoints within that region into a single daily rainfall distribution for each summer (Methods).  While the cutoff for Southeast India is negatively correlated with NINO3.4 and marginally statistically significant (\({r=-0.16}\), \(p=0.08\)), the corresponding correlations for All-Monsoonal India, the Central Monsoon Zone, and the Western Ghats regions are all significantly positive (\({r=+0.34}\), \({{+}0.34}\), \({+}0.20\) and \({p=2\times10^{-4}}\), \(2\times10^{-4}\), 0.03, respectively; Methods).  Composites combining either all summers 1901-2020, the 36 / summers, or the 41 / summers behave similarly: except for Southeast India, the region-aggregated cutoffs are all significantly larger for the / than the / composites (Fig. <ref>b; Methods).

(a) Composite regional cutoffs (in mm) in circles.
Bars are 95% confidence intervals based on yearly cutoff values, and squares are the estimates from gamma fits (Methods).
(b) ENSO risk ratios: restricting to rainy days (\(>1\) mm), the probability that rainfall exceeds a given rain rate in the / composite divided by that in the / composite.  The horizontal axis is the percentile daily rain rate across all rainy days in all years and points for that region.  Overlaid circles show each region's climatological cutoff.
(c) Same as panel c, but with logarithmically spaced horizontal axis and wider vertical axis to better discern extremes.

These signals in extremes differ qualitatively from those in less severe rain rates, which are suppressed throughout the domain by /.  Pointwise exceedance counts—the number of days in each summer exceeding a specified local rainfall percentile—are predominantly negatively correlated with NINO3.4 for low through moderate percentiles (Extended Data Fig. <ref>a-c), becoming positive for the climatologically wet regions only around the 95th percentile and above.  The same overall behavior also holds using fixed accumulations in mm (Extended Data Fig. <ref>d-f).

This can also be seen in ENSO risk ratios for each region <cit.>: at each daily accumulation, the fraction of rainy days exceeding that value in / summers divided by the corresponding fraction for / summers (Fig. <ref>b,c; Methods).  Accumulations are inhibited under / compared to / at all percentiles in Southeast India, up to the 75th percentile in the Central Monsoon Zone, and up to the 90th percentile for both the whole domain and the Western Ghats coastal band.  For the latter two especially, the percentile above which likelihoods are enhanced under / occurs near that of the climatological cutoff (overlaid circles in Fig. <ref>b), above which probabilities increase quasi-exponentially with rain rate  (Fig. <ref>c)—both signatures of an increased cutoff <cit.>.  At the 99.99th rain-intensity percentile, for example, the increases are approximately 43%, 58%, and 70% respectively for the whole domain, Central Monsoon Zone, and Western Ghats.

Cutoffs and risk ratios estimated from gamma fits to each composite distribution yield similar results (Methods; Extended Data Fig. <ref>), as do other metrics of extreme rainfall  (Methods; Extended Data Fig. <ref>).  These include correlations of block maxima with NINO3.4, pointwise Gumbel fits to the block maxima for / / / composites, and quantile regressions of daily rainfall / JJAS NINO3.4 at the 99th and 99.9th percentiles.  Results are also similar in the Indian Meteorological Department \(1^\circ\times1^\circ\) gridded daily rainfall dataset  <cit.>, which importantly for extremes <cit.> is generated from a fixed network of rain gauges (Methods; dashed curves and unfilled elements in Fig. <ref>, Fig. <ref> and Extended Data Figs. <ref> and <ref>).

Finally, while /'s propensity to generate summer-mean Indian drought notoriously weakened from the pre-satellite (1901-1978) to satellite (1979-2020) eras <cit.>—with respective NINO3.4 correlations \({r=-0.59}\) to -0.41—the ENSO signals in extremes are more robust (Extended Data Fig. <ref>; Methods).  This holds in terms of the overall geographic patterns seen in pointwise cutoffs, linear regressions of the region-aggregated cutoffs on NINO3.4, and an index linearly combining the region-aggregated cutoffs of the Core Monsoon Zone, Western Ghats, and Southeast India regions (Methods).  Nevertheless, within peninsular India both the increase in extremes with NINO3.4 in the southwestern coastal band and the decrease in the southeast weaken nontrivially.  As such for the region-aggregate cutoffs only the Central Monsoon Zone's regression against NINO3.4 is statistically significant within the satellite era.  The extent to which these modest differences between epochs are physically forced / merely sampling-driven <cit.> remains to be seen.

§ MECHANISMS: ENSO-FORCED CHANGES IN CONVECTIVE BUOYANCY EXTREMES

We expect extreme rain where near-surface air is especially warm and moist relative to air above it that, too, is very moist.  These conditions enable the near-surface air to precipitate out copious water as it ascends deeply, despite losing buoyancy to cooling through expansion and to dilution through mixing with the drier ambient air.  These processes have been encapsulated into a scalar metric of convective buoyancy known as \(B_L\) <cit.>, which we compute for each gridpoint and summer day of 1979 to 2020 using ERA5 reanalysis data <cit.> (Methods).

(a) Spearman correlation coefficient across all JJAS days 1979-2020 at each gridpoint between daily rainfall and the convective buoyancy metric.
(b) In dotted curves, PDFs of the buoyancy metric for each region.  In solid curves, daily rainfall conditionally averaged on the buoyancy metric.
(c) For the whole monsoonal India domain, composite quantiles of rainfall conditioned on the buoyancy metric as labeled.

Daily rainfall depends strongly on local convective buoyancy, with nonlinear Spearman correlation coefficient positive at every gridpoint, up to +0.71 (Fig. <ref>a) and the expected quasi-threshold relationship <cit.> in region-aggregated averages (Fig. <ref>b) and quantiles (Fig. <ref>c) of gridpoint-wise rainfall conditioned on local buoyancy (Methods).  Rainfall is typically absent or weak until buoyancy reaches a threshold value slightly below zero, above which rainfall increases sharply, rapidly consuming buoyancy with it (Fig. <ref>b).  These conditional curves rise to higher rainfall values in the rainy southwestern coastal band and Central Monsoon Zone compared to the drier southeast India, but to first approximation they are insensitive to the state of ENSO (Fig. <ref>c).

(a,b,c) Pearson correlation coefficient of NINO3.4 vs. 99th percentile daily values of the (a) full buoyancy metric (b) just the moisture term, or (c) just the undilute instability term for each JJAS.
(d,e,f) Regional ENSO risk ratios of rainy-day buoyancy, for (d) the full metric, (e) just the moisture term, or (f) just the undilute instability term as a function of return period.

Based on correlations at each gridpoint against NINO3.4, summer-mean and moderate-percentile convective buoyancy values are, like rainfall, widely suppressed in / summers (not shown).  And also like rainfall, the right-tail values are suppressed in the southeast and northwest but enhanced over much of central-eastern India  (Fig. <ref>a for 99th percentile values, with results similar for block maxima and other high percentiles, not shown).  The individual buoyancy components of moisture and undilute instability behave similarly to one another over these parts of the domain  (Fig. <ref>b,c): both very moist and very gravitationally unstable days are predominantly enhanced over eastern-central India / suppressed over the northwest and southeast.

What, then, accounts for these convective buoyancy signals?  Preliminary results point toward synoptic low-pressure systems <cit.>, which generate heavy rainfall within their cores <cit.> and track predominantly over Central India, though occasionally traveling farther south or extending to Rajasthan in the northwest <cit.>.  A prior study indicates that the core of monsoon low-pressure systems are very humid and that during / summers the rainfall within those cores is enhanced <cit.>.  And preliminary findings (Yujia You, pers. comm.) indicate that those low-pressure system tracks over and ending in central India are favored by /, at the expense of those over peninsular India and those long tracks reaching northwestern India.  If true, this would promote very humid, buoyant days in central-eastern India and suppress them in the southeast and northwest.

For the southwestern coastal band, which is less directly influenced by these low pressure systems, nearly saturated days above the boundary layer are strongly suppressed in / summers while highly unstable days are strongly enhanced.  This undilute buoyancy signal emerges to some degree for the other regions as well, as is particularly evident in region-aggregated ENSO risk ratios of the convective buoyancy terms on rainy days (Fig. <ref>d-f; Methods).  Moreover, based on differences in the ENSO composites of rainfall conditioned on the buoyancy terms (not shown), tropospheric dryness inhibits heavy rainfall somewhat less during / summers; this holds for each region but especially the southwest.  We therefore speculate that the anomalous free-tropospheric subsidence and dryness under / throughout the domain (Extended Data Fig. <ref>) suppress moderate rain events, enabling higher undilute instability values to occasionally build up, fueling extreme convection once a passing wave or advective event re-moistens the lower free troposphere just enough.  This mechanism would be broadly consistent with observational and modeling evidence for less frequent but more intense rain events in drier ambient conditions <cit.>.  Moreover, though it would in principle operate throughout the domain, plausibly it would be strongest in the southwest, where the prevailing monsoon southwesterlies reliably converging warm air into the boundary layer, promoting buoyancy.

§ OUTLOOK: SEASONAL FORECASTS, CLIMATE CHANGE, AND OTHER REGIONS

Despite ENSO's known “spring predictability barrier” <cit.>, existing seasonal forecasts of ENSO during boreal summer are skillful <cit.>; for example, the multi-model-mean forecasted JJAS value of NINO3.4 from the North American Multi-Model Ensemble (NMME) <cit.> seasonal forecast integrations initialized each May 1st are correlated with the actual value at \(r=0.86\) over 1981-2020 (Methods).  We have re-computed the correlations between the JJAS cutoff rain rate and NINO3.4 using this forecasted value, and the overall signals remain (Extended Data Fig. <ref>).  This preliminary result raises the prospect of skillful ENSO-based seasonal forecasts of extreme rainfall probability within monsoonal India.

An observed, long-term trend toward more extreme daily rainfall events in central India into the 2000s has been argued to stem from mean moistening of the troposphere due to global warming <cit.>, as have future projections of increased extreme rainfall events globally <cit.>.  Notwithstanding that more recent trends in the Indian summer monsoon have been flat or negative <cit.>, underlying this argument are the assumptions that right-tail daily moisture events (a) track with the season-mean value and (b) predominate over changes in undilute instability.  Based on our results, however, for ENSO (a) is violated within central India and (b) is violated within the southwestern coastal band—precisely the regions where daily extreme rainfall is projected to increase most in model simulations under future global warming <cit.>.  This argues for revisiting both observed trends and model projections of extreme rainfall, incorporating the full daily distributions of rainfall and of the convective buoyancy terms.  Only then can we meaningfully evaluate whether these observed, interannual, internally generated ENSO signals constitute a meaningful emergent constraint on model-projected, secular, externally forced trends <cit.>.

Our approach could also be applied to climate variability modes in the Indian <cit.> and Atlantic Oceans <cit.> thought to influence the summer monsoon, as well as to hourly rainfall.  Many Indian summer monsoon rain events span considerably less than 24 hours <cit.>, and a study <cit.> using hourly rainfall suggests that the average rain rate during rainy hours increases with / over much of the Tropics.  Perhaps, as for daily rain in the Indian summer monsoon, these signals reflect buoyancy-driven increased extremes.

§ METHODS

§.§ Region definitions

The four subregions of India are those used by <cit.>.  The “All Monsoonal India” (AMI) region comprises all gridpoints in the monsoonal India domain defined by <cit.>.  The remaining three are non-overlapping subsets of AMI.  The Central Monsoon Zone (CMZ) definition also follows <cit.>, spanning from Pakistan to the Bay of Bengal in a wide band north of peninsular India.  The Western Ghats (WG) region is defined as all land points within the polygon whose vertices are (76.5\(^\circ\)E, 7.5\(^\circ\)N), (78.5\(^\circ\)E, 7.5\(^\circ\)N), (74\(^\circ\)E, 18\(^\circ\)N), (74\(^\circ\)E, 21.5\(^\circ\)N), (72.25\(^\circ\)E, 21.5\(^\circ\)N), and (72.25\(^\circ\)E, 18\(^\circ\)N).  The Southeast India (SEI) region includes all points south of 18\(^\circ\)N that are not within WG.  The WG and SEI region borders were defined by <cit.> explicitly based on JJAS-mean rainfall properties in the Indian Meteorological Department 0.25\(^\circ\) dataset to distinguish the coherent regimes within peninsular India, but neither they nor the other regions were defined explicitly with reference to extreme rainfall characteristics.

Region-aggregated distributions are computed by combining fields from all gridpoints of a given region into a single distribution without any area weighting.  Region-averaged fields are area-weighted averages across all gridpoints within each region.

§.§ Indian Meteorological Department \({0.25^\circ\times0.25^\circ}\) rainfall dataset

Daily rainfall data comes from the 0.25\(^\circ\times\)0.25\(^\circ\) latitude-longitude gridded product produced by the Indian Meteorological Department (IMD) <cit.>.  This dataset spans all days of the year from January 1, 1901 to the present, and we use all days in June, July, August, and September (JJAS) from 1901 through 2020.  This gridded product is derived from a dense network of in situ rain gauges with coverage throughout India, using a simple interpolation procedure to go from the irregular station coverage to the regular latitude-longitude grid.  The rain gauge coverage varies in time, with \(\sim\)1500 stations in 1901, peaking above 4000 stations in the early 1990s, and then decreasing appreciably to \(\sim\)2000 in the 2010s <cit.>.  These secular trends in station coverage have been argued to strongly influence inferences about long-term trends in daily rainfall extremes <cit.>.  For our purposes, however, for interannual variability it is unlikely that the fluctuations in station coverage would alias onto ENSO.

This dataset contains spurious values in a few gridpoints and times, which appear to us to stem from errant reported values from one or more stations which then propagate in space through the interpolation procedure.  In particular, in northeast India directly along the border with Bangladesh, starting in the 1970s values become insensible.  We mask these out along with the rest of northeastern and far northern India as part of our “monsoonal India” mask <cit.>.  Despite these problems, the density of rain gauges dating back to the 19th century and the stewardship of the resulting data by the Indian Meteorological Department make this a very high-quality dataset.

§.§ Cutoff rain rate

We compute the cutoff daily rain rate using the method of moments, which for a given daily rainfall distribution is simply the variance divided by the mean.  This measure is motivated by the correspondence of daily rainfall distributions to gamma distributions <cit.>.  In terms of shape parameter \(k\) and scale parameter \(\theta\), the gamma probability density function, \(f(x)\), follows \(f(x)\propto x^{k-1}e^{-x/\theta}\), from which the cutoff is easily identified as the scale parameter \(\theta\).  Moreover, for gamma distributions the variance is exactly \(k\theta^2\) and the mean is \(k\theta\).  Thus, insofar as the daily rainfall distributions approximately follow gamma distributions, the variance-to-mean ratio usefully estimates the cutoff.

The method of moments is easily calculated, but being computed over the entire distribution, it is influenced by behavior at the low end of the daily rainfall distribution as well as the extremes.  Given the strong influence of ENSO on no- and low-rainfall day frequencies, this potentially makes it a biased metric of extremes.  We have tested the sensitivity to behaviors at low rain rates by computing it restricting to rain rates over different thresholds.  This increases the mean of the remaining sample, yielding a low-biased cutoff estimate, but results in terms of correlations with NINO3.4 or ratios of ENSO composite quantities are qualitatively insensitive (not shown).

An alternative method for computing the cutoff is via a linear regression of the daily rainfall PDF in semi-log space <cit.>.  Denoting the daily rainfall PDF as \(f(P)\), a linear regression is performed of \(\log(f(P))\) on \(P\) but restricting to the right tail of the distribution (which, being quasi-exponential, is close to a straight line); denoting the slope of this regression \(r\), the cutoff estimate is then \(-1/r\).  While this avoids contamination by the low end of the distribution, we find it strongly determined by the few largest values: for the parameter choices explored, we find the regional cutoffs for each JJAS computed via the regression to all be correlated at \(r>0.9\) with the region's JJAS block maximum rain, compared to \(r\sim0.5-0.7\) with the method of moments.

As such, we elect to use the cutoff computed from the method of moments, taking it to be more representative of behaviors throughout the right tails of the daily rainfall distributions.

§.§ NINO3.4 index

As our index of ENSO, we compute the standard NINO3.4 index of sea surface temperature (SST) anomalies averaged over the central and eastern equatorial Pacific: 120\(^\circ\)–170\(^\circ\)W, 5\(^\circ\)S–5\(^\circ\)N.  The monthly SSTs come from the National Oceanic and Atmospheric Administration Extended Reconstructed SST (ERSST) dataset, version 5 <cit.>.  The SSTs are first averaged with area-weighting over the NINO3.4 region for each month.  They are then averaged over JJAS for each year.  The linear trend computed by least squares regression over 1901-2020 of this timeseries of JJAS values is then computed and subtracted off.  Finally, the time mean of this detrended timeseries is computed and subtracted off to generate an anomaly timeseries; the result is the NINO3.4 index we use.

We have also repeated key calculations using JJAS deviations of the NINO3.4 averaged SSTs from a 30-year rolling mean with no detrending, more akin (but not identical) to the NOAA Oceanic Niño Index (ONI), finding little sensitivity.  For example, the region cutoffs for each JJAS are correlated with ONI at \({r=0.29}\), 0.29, 0.19, and -0.19 for AMI, CMZ, WG, and SEI, respectively, all less positive but within 0.05 of the corresponding correlations using NINO3.4.

§.§ Detrending

The trend for each field is computed via least-squares regression in time over 1901-2020.  De-trended fields are constructed by subtracting this trend from the full field.  All correlation coefficients are computed on fields linearly detrended in this way.

§.§ Statistical significance estimates

We use a standard bootstrapping technique to assess the statistical significance of the pointwise correlations between the cutoff and NINO3.4.  Specifically, we randomly draw with replacement 120 years from the 1901-2020 period, taking the NINO3.4 value and each point's cutoff for that JJAS.  We then compute the Pearson correlation coefficient at each gridpoint for that randomly drawn sample of 120 years.  We repeat this procedure 10,000 times, and finally construct a 95% confidence interval by selecting the 2.5th and 97.5th percentiles for each gridpoint of these 10,000 correlation coefficients.  The inferred confidence intervals converge rapidly as the number of bootstrap draws exceeds 100 (not shown).  Those points whose confidence interval do not cross zero we deem statistically significant at the 95% confidence level.

For the ENSO-composite region-aggregated cutoffs, we proceed as follows.  First, we compute region-aggregated cutoff rain rates for each individual JJAS.  For both ENSO composites and for the composite over all years 1901-2020, the average across these yearly cutoffs is similar to the cutoffs computed across all years (compare filled circles to center of corresponding vertical lines in Fig. <ref>b).  We therefore use the sample sizes (36 / years, 41 / years, 120 total years) and the interannual standard deviations in the yearly cutoffs for each region and ENSO composite to construct the 95% confidence intervals (vertical lines in Fig. <ref>b).  These regional yearly cutoffs are reasonably normally distributed (not shown), and so we employ a standard $t$-test to assess the statistical significance of the correlations between the yearly cutoffs and NINO3.4.  The resulting \(p\) values are \(p=0.0002\), 0.0003, 0.03, and 0.06 for AMI, CMZ, WG, and SEI, respectively.

§.§ ENSO composites and risk ratios

For ENSO-based composites, we use a \(\pm0.3\) K threshold of the NINO3.4 index: all summers w/ JJAS NINO3.4 exceeding +0.3 K are designated /, those between -0.3 and +0.3 neutral, and those less than -0.3 /.  These -0.3 and +0.3 threshold values correspond to the 34th and 70th percentiles of the JJAS NINO3.4 distribution respectively.  Results appear insensitive to this threshold from values of 0 to roughly \(\pm1.0\).  For example, the AMI composite cutoffs for a 0 K threshold are 36.7 and 33.8 mm respectively for / and /; for a \(\pm1\) K threshold they are 36.4 and 33.5 mm, respectively, with only small fluctuations for intermediate values.  The \(\pm0.3\) K threshold was subjectively chosen as a compromise between maximizing the sample size of each composite / maximizing the difference in the large-scale atmospheric state between the composites.

To compute the ENSO risk ratio for a given field, first we compute the empirical cumulative distribution function (CDF, denoted \(F(x)\)) for that field for the / composite and separately for the / composite.  By definition of the CDF, in either composite, the empirical probability of exceeding a given value equals one minus the CDF at that value: \(P(X\geq x)=1-F(x)\).  As such, the ENSO risk ratio is computed as \((1-F_\mr{EN}(x))/(1-F_\mr{LN}(x))\), where the subscripts “EN” and “LN” refer to / and / respectively.

We compute risk ratios in this manner for daily rain accumulation on rainy days (using a standard \({>}1\) mm threshold).  We then interpolate these risk ratios, which are functions of daily rain accumulation in mm, to climatological quantiles of daily rainfall across all years 1901-2020 and all points in the given region.  Denoting a given quantile \(q\), the quantiles are converted as standard to return periods as \(1/(1-q)\).  Because the risk ratios restrict to rainy days, the return periods correspond to JJAS rainy days as well rather than JJAS calendar days.

§.§ Exceedance counts

The exceedance count is defined for each gridpoint, year, and daily rainfall percentile as the number of days in that year's JJAS and gridpoint exceeding the climatological rainfall at the given percentile across all years and all days of JJAS.  For example, if a given year and gridpoint experienced rainfall exceeding the climatological 90th percentile on three days, the exceedance count at the 90th percentile would be three.

The exceedance count can also be applied to fixed / thresholds rather than local percentiles.  Extended Data Fig. <ref>d,e,f shows the correlations / NINO3.4 of the yearly exceedance counts for three of the official thresholds in / the Indian Meteorological Department uses to categorize rainfall severity.

§.§ Explicit gamma fits

For each region and each ENSO-composite daily rainfall distribution, a gamma distribution was fitted for all daily accumulations \(>0.1\) mm using the  function from the SciPy package <cit.> for Python.  This yields for each distribution values of the three parameters of the gamma distribution, location, shape, and scale.  In the exact gamma limit, the cutoff as estimated by the method of moments is identical to the scale parameter, and we therefore take the scale parameter for each of these fits as the gamma-estimated cutoff.  The ENSO risk ratios computed for the explicit gamma fits were computed from their CDFs identically to those for the empirical distributions.

§.§ Other extreme daily rainfall metrics

The block maximum daily rainfall for each JJAS is simply the single daily maximum rainfall in that JJAS at each gridpoint.  We create / and / composite block-maxima distributions for each gridpoint by combining the block maxima across all / summers or all / summers into a single distribution.  We then fit a Gumbel distribution to each / composite and / composite using the  function, yielding for each gridpoint and composite a scale parameter value and a location parameter value.  Increasing either the scale or location parameter of the Gumbel distribution increases right-tail probabilities.

Quantile regression determines the linear slope of the dependent variable (in this case, daily rainfall) / the independent variable (in this case, the JJAS NINO3.4 value) such that, for the specified quantile \(q\), a fraction \(q\) of the dependent variable points fall below that line and a fraction \(1-q\) of points fall above that line <cit.>.  For a given JJAS, each day is assigned the JJAS-mean NINO3.4 value, and for each gridpoint all JJAS days of all years are combined into one distribution.  We perform quantile regression for the 99th and 99.9th percentiles.

§.§ Indian Meteorological Department \({1^\circ\times1^\circ}\) rainfall dataset

We have replicated all key analyses on the Indian Meteorological Department's 1\(^\circ\times\)1\(^\circ\) gridded product <cit.>.  This product uses the same interpolation procedure of Indian Meteorological Department station data, but in addition to being coarser, it differs from the 0.25\(^\circ\) product in that its network of 1,803 stations is fixed in time.  There are no substantive discrepancies in our key findings between the two products.  Perhaps most notable is that daily rainfall values within the very narrow WG coastal band reach considerably larger values in the higher-resolution product, which is reflected in the WG region-aggregated daily cutoff (Fig. <ref>b).  But this amounts to a approximately constant offset; the ENSO signals in this and other regions are similar between the two products.

§.§ Pre-satellite and satellite era results

A multiple linear regression over 1901-2020 on NINO3.4 of the Central Monsoon Zone, Western Ghats, and Southeast India cutoffs yields respective weights of +1.0, +0.4, and -1.0, and the resulting timeseries is correlated with NINO3.4 for the full, pre-satellite, and satellite periods respectively at \({r=0.45}\), 0.46, and 0.44.  In the satellite era, the JJAS NINO3.4 \(\pm\)0.3 K threshold results in 14 / summers and 13 / summers.

§.§ ERA5 reanalysis data

Large-scale environmental fields come from the European Center for Medium-range Weather Forecasting (ECMWF) ERA5 reanalysis dataset <cit.>, which provides data at hourly and monthly resolutions on a global 0.25\(\times\)0.25\(^\circ\) grid.  Conveniently, over India this grid aligns identically with that of the Indian Meteorological Department rainfall.  Though the ERA5 dataset now extends from 1950 to present, we restrict to the satellite era, 1979-2020, as the amount of observational data assimilated into the underlying numerical model increases dramatically with the availability of satellite retrievals.  We compute JJAS-mean quantities from ERA5 monthly fields for June through September, and we compute daily-mean fields from ERA5 hourly fields, averaging over all 24 hourly fields.

§.§ \(B_L\) convective buoyancy metric

We compute the convective buoyancy metric, \(B_L\), following <cit.> as

\begin{equation}
\label{eq:bl-def}
B_L\equiv g\left[\wpbl\frac{\thetaepbl-\thetaesatlft}{\thetaesatlft}+\wlft\frac{\thetaelft-\thetaesatlft}{\thetaesatlft}\right],
\end{equation}

where \(g\) is gravity, \(\theta_e\) is equivalent potential temperature, \(\thetaepbl\) is \(\theta_e\) averaged over the boundary layer, \(\theta_e^*\) is saturation equivalent potential temperature, \(\thetaelft\) is \(\theta_e\) averaged over the lower free troposphere, \(\thetaesatlft\) is \(\theta^*_e\) averaged over the lower free troposphere, and \(\wpbl\) and \(\wlft\) are the weights given to the undilute instability and subsaturation terms respectively.  Equivalent potential temperature is computed conventionally as

\begin{equation}
\label{eq:theta-e}
\theta_e\equiv T\left(\frac{p_0}{p_\mr{d}}\right)^{R_\mr{d}/c_p}\mathcal{H}^{-R_\mr{v}r_\mr{v}/c_p}\exp\left(\frac{L_\mr{v}r_\mr{v}}{c_p T}\right),
\end{equation}

where \(T\) is temperature, \({p_0=1000}\) hPa is a reference pressure, \(p_\mr{d}\) is the dry air partial pressure, \(R_\mr{d}\) is the dry air gas constant, \(c_p\) is the specific heat of dry air at constant pressure, \(\mathcal{H}\) is relative humidity, \(R_\mr{v}\) is the gas constant of water vapor, and \(r_\mr{v}\) is the water vapor mixing ratio.  This definition neglects the contribution to \(\theta_e\) of liquid water.  Saturation equivalent potential temperature is computed using (<ref>) but setting \({\mathcal{H}=1}\) and the vapor mixing ratio to its saturation value.

The layer weights are given by

\begin{align}
\label{eq:wb-def}
\wpbl&\equiv \frac{a\Delta p_\pbl}{b\Delta p_\lft}\ln\left(\frac{a\Delta p_\pbl + b\Delta p_\lft}{a\Delta p_\pbl}\right)\\\
\wlft&\equiv 1 - \wpbl,
\end{align}

where \(a\) and \(b\) are parameters relating to the relative mass inflow rates in the PBL and LFT layers, respectively, and \(\Delta p_\pbl\) and \(\Delta p_\lft\) are the pressure thicknesses of the PBL and LFT layers, respectively.  We follow <cit.> in setting \(a=b=0.5\), but we compute \(\Delta p_\pbl\) and \(\Delta p_\lft\) for each day and gridpoint as follows.

Following convention <cit.>, we define the boundary layer to extend from the local surface pressure to 150 hPa above.
To better resolve the layer boundaries and thus depths, ERA5 fields are linearly interpolated in pressure from their native resolution (of 25 hPa spacing from 1000 to 750 hPa, 50 hPa spacing from 750 to 250 hPa, and 25 hPa spacing from 250 to 100 hPa) to 5 hPa spacing up to 200 hPa.  The LFT is taken to span from the next level vertically above to the last level below the freezing level, defined as where the temperature drops below 273 K.

When computing the free-tropospheric subsaturation term, \(\blsubsat\), and the subcloud undilute instability term, \(\blinstab\), individually, the layer-depth weights are included in each:

\begin{align}
\label{eq:bl-terms}
\blsubsat&\equiv g\wlft\frac{\thetaepbl-\thetaesatlft}{\thetaesatlft}\\
\blinstab&\equiv g\wpbl\frac{\thetaelft-\thetaesatlft}{\thetaesatlft}.
\end{align}

Given our focus on extremes, we use this version of \(B_L\) in terms of equivalent potential temperatures rather than the simplified version in terms of moist enthalpies <cit.>, as the latter loses some of the variations within either tail.  Our use of a spatiotemporally varying PBL bounds and LFT layer thickness defined in terms of the local surface pressure and freezing level are also more involved than the original version which uses constant values of 1000 hPa as the surface pressure, 850 hPa as the PBL top, and 500 hPa as the LFT top.

§.§ ENSO risk ratios of \(B_L\)

Risk ratios for \(B_L\), \(\blsubsat\), and \(\blinstab\) components are computed identically to those for rainfall, other than being restricted to the satellite era 1979-2020.  This includes restricting to rainy days and the risk ratios then being interpolated to regional climatological quantiles and then return periods.

§.§ Data availability statement

Indian Meteorological Department rainfall data used is available at <https://www.imdpune.gov.in/cmpg/Griddata/Rainfall_25_NetCDF.html> and <https://www.imdpune.gov.in/cmpg/Griddata/Rainfall_1_NetCDF.html>.  ERSST data is available at <https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/netcdf/>.  ERA5 data is available at <https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5>.  NMME data is available at <https://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/>.

S.A.H. acknowledges funding from NSF award AGS-2123327.  S.A.H., M.B., M.A.C., and A.H.S. acknowledge support from the Monsoon Mission Project under India’s Ministry of Earth Sciences.

Author contributions
S.A.H. designed the study, oversaw preliminary analyses performed by D.M.Z., performed most analyses, and led writing the manuscript.  D.M.Z. performed the initial analyses establishing the ENSO-extreme rainfall signals.  A.H.S., M.B., and M.A.C. provided scientific feedback throughout the project.   M.K.T. provided guidance on statistical procedures.  F.A. provided guidance on convective buoyancy.  All authors contributed to revising the manuscript.

Competing interests
The authors declare no competing interests.

Materials & Correspondence
All correspondence and material requests should be addressed to S.A.H.

§ EXTENDED DATA FIGURES

In color shading according to each colorbar, (a) climatological JJAS-mean rainfall (units mm/day), (b) JJAS-mean rainfall interannual standard deviation (units mm/day), (c) daily rainfall standard deviation (units mm), (d) average number of days exceeding the Indian Meteorological Department's “heavy rain” threshold of 64.5 mm, (e) average across years of each year's JJAS block maximum rain (units mm), and (f) the climatological 99.9th percentile daily accumulation across all JJAS days (units mm).  For each, brown shades denote values less than approximately the median across gridpoints, and blue shades denote values above that value.

Timeseries of daily rainfall block maximum within each JJAS at three representative points.  In very rainy points of WG (dark green), there are very heavy rain events nearly every summer; in dry points of SEI (brown), the heaviest events are much weaker.  In Gujarat (far western CMZ; orange), most summers see modest max rainfall but occasionally receive extremely high values.

Scatterplot of the JJAS cutoff-NINO3.4 correlation coefficient / climatological JJAS-mean rainfall for points with JJAS climatology \(<12\) / restricting to points below the JJAS climatological rain 95th percentile across gridpoints (12.0 /).  Overlaid lines are the corresponding linear regressions, from the \(0.25\times 0.25^\circ\) dataset in solid and the \(1^\circ\times 1^\circ\) dataset in dashed.

In color shading according to the colorbar, correlation coefficient of exceedance counts for each JJAS against NINO3.4 of the local climatological (a) 95th, (b) 99th, and (c) 99.9th percentile, or of the Indian Meteorological Department's thresholds for (d) moderate, (e) heavy, or (f) very heavy daily accumulations, with the corresponding value in mm printed in each label.

(a) In squares, composite regional cutoffs computed from gamma fits; all other plotted elements are identical to Fig. <ref>(a).  (b) ENSO risk ratios of daily rainfall for each of the four monsoonal India regions according to the colors and labels.  Dotted curves are the actual risk ratios from the $0.25^\circ\times0.25^\circ$ Indian Meteorological Department rainfall dataset.  Solid curves are the corresponding gamma fits.  Dashed curves are the gamma fits for the $1^\circ\times1^\circ$ dataset.

(a) Pearson correlation coefficient between the JJAS block maximum daily rainfall and NINO3.4, 1901-2020.  (b and c) Ratio of / to / composite Gumbel fits to JJAS block maximum daily rainfall at each gridpoint, for the (b) location and (c) scale parameters.  (d and e) Quantile regression for the (c) 99th and (d) 99.9th percentile daily rainfall (including both rainy and dry days) against NINO3.4.

(a and b) Pointwise correlations between the JJAS cutoff and NINO3.4 restricting to (a) the pre-satellite era of 1901-1978 and (b) the satellite era of 1979-2020.  (c) ENSO risk ratios for each region restricting to (dashed) the pre-satellite era and (solid) the satellite era.  (d through g) Scatterplots of the region-aggregated cutoff in each JJAS / NINO3.4, for (d) All Monsoonal India, (e) Central Monsoon Zone, (f) Western Ghats, and (g) Southeast India.  Overlaid lines are corresponding linear regressions.  Yellow elements are for the pre-satellite era, dark red elements are for the satellite era, and black elements are for the full period.  Printed values above each panel are the regression slopes for that region, in mm K$^{-1}$, with values in bold statistically significant at the \(p=0.05\) level.  Light red and blue shading denote / and / conditions based on the \(\pm0.3\) K threshold used elsewhere for constructing composites.  Note differing vertical axis spans across d-f.

JJAS-mean ENSO composite anomaly profiles, / minus /, averaged over the AMI domain, in relative humidity and vertical velocity.

For the 1980-2020 period for which NMME data is available, correlation coefficient of the local JJAS cutoff / the JJAS NINO3.4 value computed (a) directly from ERSST, or (b) forecasted by the NMME multi-model mean of runs initialized at the preceding May 1st.  Panel (c) shows their difference.


# Infinite-Dimensional Adaptive Boundary Observer for Inner-Domain Temperature
Estimation of 3D Electrosurgical Processes using Surface Thermography Sensing

Hamza El-Kebir1, Junren Ran2, Martin Ostoja-Starzewski3, Richard Berlin4,
Joseph Bentsman2, and Leonardo P. Chamorro2 1H. El-Kebir is with the
Department of Aerospace Engineering, University of Illinois Urbana-Champaign,
Urbana, IL, 61801 USA<EMAIL_ADDRESS>Ran, J. Bentsman, and L. P.
Chamorro are with the Department of Mechanical Science and Engineering,
University of Illinois Urbana-Champaign, Urbana, IL, 61801 USA
<EMAIL_ADDRESS>Ostoja-Starzewski is with the Department of
Mechanical Science and Engineering, the Beckman Institute, and the Institute
for Condensed Matter Theory, University of Illinois Urbana-Champaign, Urbana,
IL, 61801 USA<EMAIL_ADDRESS>Berlin is with the Department of Trauma
Surgery, Carle Hospital and the University of Illinois Urbana-Champaign,
Urbana, IL 61801 USA.Research reported in this publication was supported by
the National Institute of Biomedical Imaging and Bioengineering of the
National Institutes of Health under award number R01EB029766. The content is
solely the responsibility of the authors and does not necessarily represent
the official views of the National Institutes of Health.

###### Abstract

We present a novel 3D adaptive observer framework for use in the determination
of subsurface organic tissue temperatures in electrosurgery. The observer
structure leverages pointwise 2D surface temperature readings obtained from a
real-time infrared thermographer for both parameter estimation and temperature
field observation. We introduce a novel approach to decoupled parameter
adaptation and estimation, wherein the parameter estimation can run in real-
time, while the observer loop runs on a slower time scale. To achieve this, we
introduce a novel parameter estimation method known as attention-based noise-
robust averaging, in which surface thermography time series are used to
directly estimate the tissue’s diffusivity. Our observer contains a real-time
parameter adaptation component based on this diffusivity adaptation law, as
well as a Luenberger-type corrector based on the sensed surface temperature.
In this work, we also present a novel model structure adapted to the setting
of robotic surgery, wherein we model the electrosurgical heat distribution as
a compactly supported magnitude- and velocity-controlled heat source involving
a new nonlinear input mapping. We demonstrate satisfactory performance of the
adaptive observer in simulation, using real-life experimental ex vivo porcine
tissue data.

##  I Introduction

Over the past decades, the advantages of robot-assisted surgery over
conventional laparoscopy have become increasingly apparent [1]. As the robotic
surgery performance demands grow, so does the need for a deeper understanding
of the physical phenomena governing both the tools and the tissue treated, as
well as the control laws for the precise attainment of the surgical
objectives.

Electrosurgery relies on the use of high power density radio frequency
currents to actively heat organic tissue, allowing it to be denatured,
coagulated, desiccated, fulgurated, or incised [2]. One of the key advantages
of electrosurgery is its capability of simultaneous cutting and coagulation,
providing blood stoppage for complex surgical tasks. Since this technique
allows for precise ablation of tissue with very little collateral damage
($\sim$100–400 $\mu$m), it is commonly used in practice, with over half of the
surgical procedures employing it [2]. An illustration of the electrosurgical
process along with the observer setup considered in this work is shown in Fig.
1.

Figure 1: Illustration of the electrosurgical process as it relates to the
adaptive inner-domain temperature observer presented in this work. Given the
electrosurgical probe power and position commands, as well as a surface
thermography image, we estimate the thermal diffusivity and employ an adaptive
observer law to obtain the inner-domain temperature field.

Electrosurgical control can be posed as a boundary control problem
characterized by a controllable heat flux, as dictated by the power setting,
and cathode position actuation. In previous work [3, 4], the authors have
considered a 1D setting in which the tissue thermodynamics are modeled using a
Stefan problem formulation, where a homogeneous substance undergoing a moving
phase change from virgin to a denatured tissue is considered. In [3] a
stationary probe was considered, with a novel approach to non-collocated
thermal sensing being introduced. The model in [3] was applied as part of a
real-life planar cutting task using a hypermodel-based model predictive
control algorithm in [4]. However, to guarantee a tight attainment of the
control objective in the case of a moving probe location, which is essential
in practice [2], considering the process in 3D becomes a necessity. In this
work, we consider a spatially distributed controlled moving heat source, while
accounting for chemical reactions by incorporating a bounded disturbance. In
this setting, as the first approximation, a homogeneous Fourier heat equation
is employed, with a moving heat volumetric source that has compact spatial
support. In this work, we develop a 3D observer structure based on surface
temperatures with the goal of estimating the subsurface temperature
distribution; the latter goal is essential in assessing and preventing
subsurface tissue damage. In practice, these inner-domains temperatures cannot
be sensed in a noninvasive fashion, thereby significantly limiting any model-
based control efforts for safe autonomous electrosurgery. In addition to
proving convergence under bounded heat input disturbances in the case of known
thermophysical parameters, we also introduce a novel real-time parameter
adaptation loop in which _only_ surface temperature measurements are used.
This provides the first adaptive surface-temperature-based observer for
electrosurgical processes.

In the context of biomechanics and electrosurgery, the first works dealing
with temperature distributions in nonhomogeneous tissue models considered the
thermal response of layered media when subjected to microwave excitation using
microwave radiometry [5], building on earlier work that considered homogeneous
tissue models [6]. In all of these works, experimental verification rested on
the use of surface temperatures measurements; indeed, obtaining inner-domain
temperatures is not possible without invasively altering the heat transfer
process, e.g., by inserting thermocouples. Recent work [7] has used multi-
angle surface thermography images to estimate the tomographic, or volumetric,
temperature field. In [7], an application-specific deep neural network model
was trained to deduce the inner-domain temperatures from the surface
temperatures after being trained on synthetic data based on an unforced
Fourier heat transfer model with constant thermal diffusivity throughout the
entire domain. Given this data, such a model is only capable of estimating
inner-domain temperatures for experiments that involve a free thermal response
and have the same thermophysical parameters, boundary conditions, and
geometry. In the context of electrosurgery, it is essential that an external
heat input be accounted for, as well as to have an ability to obtain
subsurface temperatures for different tissue types and geometries. It is then
natural to consider the infinite-dimensional observers based on the surface
temperature feedback, as is done in this work.

In constructing model-based observers, it is essential to have a sufficiently
accurate model of the thermodynamics. While the thermal properties of the
tissue can readily be obtained in a laboratory environment using differential
scanning calorimetry (DSC) and thermogravimetric analysis (TGA) [3], these
methods are costly and time-consuming, and only provide localized estimates of
the tissue properties. More importantly, such methods are necessarily
invasive, since they require biopsic samples to be extracted which can often
only be used for one form of analysis at a time. Circumventing these issues,
which extend to nondestructive testing of materials as well, in the context of
subsurface material characterization a technique known as dynamic thermal
tomography (DTT) has been introduced by Vavilov _et al._ [8], and has recently
seen use in medical applications [9, 10]. In DTT, thermographic images are
captured at different areas on the surface, with the goal of estimating local
thermophysical characteristics (thermal conductivity, heat capacity, heat
transfer coefficients, etc.) by exciting dynamic thermal responses using
embedded collocated heating elements. Using the obtained volumetrically
varying parameter estimates, it is possible to formulate a model PDE and
obtain the subsurface temperatures as was done in [9]. An early approach to
this method for use in biomedical applications is presented in [11], but it
relies on Fourier series expansions. A major drawback of DTT is the fact that
it is an invasive tomography method, unlike other tomographic imaging
techniques such as magnetic resonance imaging. In this work, we introduce a
non-invasive real-time parameter adaptation framework based solely on the
surface temperature. To enable the latter we propose a novel thermography-
based parameter estimation method which we call _attention-based noise-robust
averaging_ (ANRA).

The paper has the following structure. The boundary control problem of
interest along with the related mathematical model are presented in Sec. II.
Sec. III presents a non-collocated surface-temperature-based observer
structure for electrosurgical action that is asymptotically convergent in the
case of known thermophysical parameters. Sec. IV presents a novel adaptive
observer structure, in which the thermophysical parameters are directly
estimated from the real-time surface thermography. Sec. V compares the
performance of both observer structures through numerical simulation.
Conclusions are drawn in Sec. VI. Formal proofs will be presented in a
separate publication.

##  II Problem Formulation

In this work, we consider the following thermodynamic model on an
$n$-dimensional compact domain $\Omega\subseteq\mathbb{R}^{n}$, based on a
standard model from the literature [12, Eq. 1]:

$\dot{z}(t)=\begin{bmatrix}\dot{x}(t)\\\
\dot{p}(t)\end{bmatrix}=\begin{bmatrix}a\nabla^{2}x(t)+T_{p(t)}u(t)q+w(t)\\\
v(t)\end{bmatrix},$ (1)

where $a\in\mathbb{R}_{+}$ denotes the thermal diffusivity, and $\nabla^{2}$
is the Laplacian operator. $x(t)$ denotes the temperature field at time $t$,
such that $x(t)\in H^{1}(\Omega)$; by $x(t,\eta)$ we denote $x(t)(\eta)$ for
$\eta\in\Omega$. The heat source position is denoted by $p(t)\in\Omega$, which
is used as a parameter in the functional $T_{p}$ defined as

$\begin{split}T_{p}:H^{1}(\Omega)&\to H^{1}(\Omega),\quad\forall
p\in\Omega,\\\
T_{p}q(\eta)&\mapsto\chi_{\Omega}(\eta)q(\eta-p),\quad\forall\eta\in\Omega,q\in
H^{1}(\Omega),\end{split}$ (2)

where $\chi_{\Omega}$ is the indicator function on the set $\Omega$. In our
model, $q\in H^{1}(\Omega)$ denotes a normalized heat source term, which is
scaled by the control input $u(t)\in\mathbb{R}_{+}$. The heat input is
perturbed by $w(t)\in\mathcal{W}:=\\{w\in
H^{1}(\Omega):\|w\|_{L^{2}(\Omega)}\leq\epsilon_{w,2},w\in
C^{0,\gamma}(\Omega)\\}$, which is bounded in the $L^{2}(\Omega)$-norm by
known $\epsilon_{w,2}\geq 0$, and is $\gamma$-Hölder continuous for
$\gamma=1/2$. Finally, $v\in C^{0}(\mathbb{R}^{n})$ is the known controlled
heat source velocity, which is such that $p(t):=p(0)+\int_{0}^{t}v(\tau)\
\mathrm{d}\tau\in\Omega$ for all $t\in[0,\infty)$. We denote the full state by
$z\in H^{1}(\Omega)\times C^{1}(\Omega)$.

In this work, we assume having access to thermographer data captured on one of
the boundaries of $\Omega$, which we denote by
$\Gamma:=\\{\eta\in\Omega:\eta_{n}=0\\}$. In a three-dimensional setting, this
would be the top surface of a slab of material. Furthermore, we consider that
temperatures are sampled on a finite Cartesian grid with uniform spacing
$\Delta\eta>0$, denoted by $\hat{\Gamma}$. In this work, we consider
$\mathrm{dim}\ \hat{\Gamma}=M<\infty$, such that
$\hat{\Gamma}:=\\{\eta^{(i)}\\}_{i=1}^{M}$.

Let $C:H^{1}(\Omega)\to\mathbb{R}^{M}$ be the output map, defined as
$Cx(t):=\left(x(t,\eta^{(i)})\right)_{i=1}^{M}$. The output is denoted by
$y(t):=Cx(t)$. We introduce a map $C^{*}$ from $\mathbb{R}^{M}$ to
$H^{1}(\Omega)$, which will be used as part of a Luenberger observer structure
in the following section:

$\left(C^{*}y(t)\right)(\eta):=\sum_{i=1}^{M}y_{i}\delta(\eta-\eta^{(i)})=\sum_{i=1}^{M}x(t,\eta^{(i)})\delta(\eta-\eta^{(i)}),$

where $\delta$ is the Dirac delta. In the remainder of this work, we make the
following assumptions:

  1. (A1)

The initial condition $x(0)\in H^{1}(\Omega)$ has known constant temperature,
i.e., $x(0,\eta)=x_{0}\chi_{\Omega}(\eta)$ for all $\eta\in\Omega$ for some
known $x_{0}\in\mathbb{R}$.

  2. (A2)

$\inf u(t)\geq 0$ and $\sup u(t)<\infty$.

  3. (A3)

On the boundary of $\Omega$, $\partial\Omega$, a trivial Neumann boundary
condition holds, i.e., $\nabla x(t)|_{\partial\Omega}\equiv 0$.

  4. (A4)

The disturbance $w(t)$ is bounded in the $L^{2}$-norm
$\|w(t)\|_{L^{2}(\Omega)}\leq\epsilon_{w,2}$ for all $t\in[0,\infty)$, for
known $\epsilon_{w,2}\geq 0$. In addition, $w(t)$ is $\gamma$-Hölder for
$\gamma=1/2$, with the Hölder coefficient being upper-bounded by a known
constant $C_{w}\geq 0$. Finally, $w\in C^{0}(\mathcal{W})$, i.e., $w$ is
continuous in time.

  5. (A5)

Heat source velocity $v\in C^{0}(\mathbb{R}^{n})$ is continuous and bounded,
such that $p(t):=p(0)+\int_{0}^{t}v(\tau)\ \mathrm{d}\tau\in\Omega$ for all
$t\in[0,\infty)$.

As can be seen in (1), we assume that there are no errors in the probe
position. In the remainder of this work, we consider control input pair
$(u(t),v(t))$ to be given for all $t\in[0,\infty)$; we also assume that there
is no error in the control inputs.

######  Remark 1.

It will be shown in Sec. V that Assumption (A1) does not necessarily need to
hold. Assumption (A2) requires boundedness and non-negativity of the control
signal, which follows from physical limitations in electrosurgery (i.e., heat
can only be supplied by the electrosurgical pencil, and only a finite power
can be supplied). Assumption (A3) is standard and can be argued to hold for
sufficiently large domains, where the far-field temperature is constant (i.e.,
in thermodynamic equilibrium).

Assumption (A4) puts constraints on the heat source disturbance. First, it
should be noted that the $L^{2}$-norm bound imposed by $\epsilon_{w,2}\geq 0$
requires that the disturbance signal is of bounded energy in space, but not in
time. This implies that the disturbance could be persistent across the entire
system run, with no ultimate boundedness being imposed. Furthermore, the
$1/2$-Hölder coefficient bound $C_{w}\geq 0$ imposes a constraint on the
sharpness of the disturbance; in particular, this constraint is more general
than a Lipschitz bound. While white noise cannot be captured as part of this
assumption, wide-band non-decaying signals can be. This $L^{2}$-norm bound can
be found by considering unmodeled dynamics such as tissue perfusion, and
accounting for the magnitude of these terms based on [12].

##  III Observer Design

We consider the following observer system:

$\dot{\hat{x}}(t)=a_{0}\nabla^{2}\hat{x}(t)+T_{p(t)}u(t)q+C^{*}G(t)(y(t)-\hat{y}(t)),$
(3)

where $G(t)\in\mathbb{R}^{M\times M}$ is the time-varying observer gain, which
is diagonal and positive definite for all $t\geq 0$, such that at least one of
the $g_{i}(t):=G_{ii}(t)$ satisfies

$\begin{split}g_{i}(t)\geq\ &\epsilon_{w,2}\ m(\Omega)^{1/2}\
\tilde{x}_{s,i}^{-2}(t)\\\
&\times\left(\max_{\eta\in\Omega}\min_{i=1,\ldots,M}\left[|\tilde{x}_{s,i}(t)|+C_{\gamma}(t)\|\eta^{(i)}-\eta\|^{\gamma}\right]\right),\end{split}$
(4)

where $m(\Omega)$ is the volume of the domain $\Omega$, and

$C_{\gamma}(t)=\frac{C_{w}}{na}(\exp(na\ t)-1).$ (5)

Note that $a_{0}$ is assumed to be known; the case of unknown thermal
diffusivity, is treated in the next section. Here, the disturbance Hölder
coefficient bound $C_{w}$ and the $L^{2}$-bound $\epsilon_{w,2}$ are known
based on Assumption (A4). Note that, as mentioned in the previous section, we
assume that $p\in C^{1}(\Omega)$ is known without error. We start by
considering asymptotic error convergence for the case where $a_{0}=a$.

### Observer Error Convergence

######  Theorem 1.

Consider the observer given by (3)–(5). Provided that assumptions (A1)–(A4)
hold, the temperature estimation error converges to zero.

######  Remark 2.

The lower bound on the gain given by (4) and (5) can be computed directly
given the internal observer state $\hat{x}(t)$, the observed surface
temperature $y(t)$, and the known bounds $\epsilon_{w,2}$ and $C_{w}$ from
(A4). Given the assumptions, $\tilde{x}_{s,i}(t)$ is continuous, which in turn
implies that the lower bound on the gain, (4), is continuous.

Given the evolution of $C_{\gamma}(t)$ from (5), it can be seen that the lower
gain bound is increasing in time, but is bounded for finite $t$. In practice,
a reasonably large initial choice of gain will result in (4) being satisfied
during some time interval $[0,T]$ for $T\in\mathbb{R}_{+}$.

Theoretically, the elements $g_{i}(t)$ could be taken arbitrarily large, since
no sensor noise is considered as part of this work; this would result in an
instantaneous correction of the temperature error at each sensing location. In
practice, one would not want to increase the gain too much, since this results
in amplification of unwanted sensing noise.

##  IV Adaptive Observer Design

We now consider the case where $a_{0}\neq a$, with unknown $a$. Given that the
true diffusivity is unknown, we introduce the following adaptive observer
structure:

$\begin{split}\dot{\hat{x}}(t)&=\bar{a}(t)\nabla^{2}\hat{x}(t)+T_{p(t)}u(t)q+C^{*}G(t)(y(t)-\hat{y}(t)),\end{split}$
(6)

where $G(t)$ obeys the conditions of Thm. 1, and $\bar{a}(t)$ is defined as

$\bar{a}(t):=a_{0}+L\int_{0}^{t}(\hat{a}(\tau)-\bar{a}(\tau))\text{d}\tau,$
(7)

where $\hat{a}(t)$ is the instantaneous thermal diffusivity estimate obtained
using the _attention-based noise-robust averaging_ (ANRA), introduced next.
The parameter adaptation gain is denoted by $L\geq 0$. In this formulation,
the parameter estimate $\bar{a}(t)$ in effect acts as an integral correction
term to the Laplacian term.

Figure 2: Overall _attention-based noise robust averaging_ (ANRA) model
architecture: (I) real-time thermographic imagery is acquired; (II) gradients
that appear in the PDE model are computed using noise-robust gradient
operators; (III) computed gradients are used to directly estimate parameters;
(IV) an attention layer based on the rate of temperature change (RTC) is used
to compute weighted (spatially averaged) parameter estimates; (V) based on a
confidence threshold on the RTC-field, final parameters are computed. Figure
3: Block diagram of the proposed adaptive observer (6).

### Adaptive Noise-Robust Averaging

To obtain $\hat{a}(t)$, we note that we may only use the surface temperature
readings from the thermographer. To use the latter for our purpose we propose
a novel framework further referred to as _attention-based noise-robust
averaging_ (ANRA). ANRA relies on noise-robust gradient operators that allow
for computing the Laplacian and time derivative while suppressing high-
frequency noise, which would ordinarily corrupt any further data processing
steps. We denote the thermal diffusivity estimate obtained by ANRA based on
discrete-time surface temperature measurements by $\hat{a}[k]$, where $k$ is
the current time sampling index.

Since our parameter estimation architecture depends on direct computation of
the spatial and temporal temperature gradients, we wish to suppress the
effects of high frequency noise, which is observed both spatially and
temporally in thermographic data. For this reason, a _noise-robust filter_ ,
i.e., one that subdues high frequency noise while leaving low frequency
content unattenuated, is desirable [13]. Pavel Holoborodko introduced noise-
robust gradient operators in backward differencing form [13], as well as
linear filters for spatial gradients [14]; this family of operators is
instrumental in our work. We denote noise robust finite difference operators
by $\hat{\partial}$.

As part of the ANRA framework, we introduce the _rate of thermal change_ (RTC)
field indicating the areas of the temperature field exhibiting the most
activity. The RTC layer serves to indicate the points likely to be noise-
driven, and allowing us to favor the high thermal activity regions. For this
reason, we propose to draw the attention of the parameter estimation algorithm
to areas of high RTC, using a data-driven _attention layer_. The RTC field is
defined as follows:

$\begin{split}\text{RTC}_{i,j}[k]&\leftarrow\exp(\beta|\hat{\partial}_{t}\bar{x}_{i,j}[k]|)-\min_{i,j}\exp(\beta|\hat{\partial}_{t}\bar{x}_{i,j}[k]|)\\\
\text{RTC}[k]&\leftarrow\text{RTC}[k]/\sum_{i,j}\text{RTC}_{i,j}[k],\end{split}$
(8)

where $k$ denotes the sample index in time, $\bar{x}[k]$ denotes the surface
temperature matrix captured by the thermographer, and $\beta\geq 1$ is a
tuning parameter that increases the RTC value separation between areas of low
and high rate of change as it increases; in this work, $\beta=100$. The choice
for this definition is as follows. We take the absolute value of the
temperature gradient at each location and exponentiate it, so as to assign a
larger positive value to a large thermal rate of change, this value scaling
nonlinearly in the gradient. We then subtract the smallest element to get the
minimum to be zero. Finally, we divide by the sum of all elements, to obtain a
matrix whose elements sum to 1. This latter property is instrumental as shown
next.

Having obtained the RTC field, we may now introduce a weighted averaging step.
The instantaneous thermal diffusivity estimate is obtained as follows:

$\hat{a}[k]=\sum_{i,j}\text{RTC}_{i,j}[k]\left\\{\hat{\partial}_{t}\bar{x}[k]/[(\hat{\partial}_{11}+\hat{\partial}_{22})\bar{x}[k]]\right\\}_{i,j}.$
(9)

In this step, we obtain a weighted average that is governed by the magnitude
of the rate of temperature change. This is motivated by the observation that
the sensor data from the areas that experience little thermal change are often
dominated by noise. In light of (1), thermal diffusivity estimates obtained
using ANRA are only valid when the power input, $u(t)$, is zero. During
surgery, such periods of inactivity frequently occur to assess the quality of
an incision that was just made, thereby providing ample opportunity to compute
$\hat{a}$ during these times. An overview of the ANRA method is shown in Fig.
2, and the full adaptive observer structure is summarized in Fig. 3.

It must be noted that when obtaining the diffusivity estimates, temperature
gradients in the depthward direction cannot be sensed. This is natural, since
thermography cannot determine temperatures _inside_ a material. This does,
however, lead to a discrepancy in the thermal diffusivity when compared to a
correct estimate based on the full depthward temperature information. We
propose a constant scaling factor on $\hat{a}[k]$ to correct for this fact, as
outlined in the following lemma.

######  Lemma 1.

Consider an infinitesimal 3D Cartesian volume of a homogeneous material, whose
temperature is governed by the heat equation $\dot{\xi}(t)=a\nabla^{2}\xi(t)$.
Let the top surface of this volume be fully insulated (i.e.,
$\partial^{2}_{\hat{n}_{\text{top}}}\xi(t)=0$). Assume that along each
Cartesian axis, the second derivative of temperature along that axis is the
same for each direction. Then, if the bottom surface is also modeled to be
insulated, while the true system is not insulated at the bottom, the effective
thermal diffusivity $a^{\prime}$ will be equal to $\frac{5}{4}a$.

Using the result from Lemma 1, it suffices to scale any thermal diffusivity
estimate obtained from surface temperature measurements by $4/5$.

######  Conjecture 1.

Given the adaptation law introduced in (6), the adaptive observer converges
asymptotically, with faster convergence compared to the non-adaptive observer
of Thm. 1.

We demonstrate Conjecture 1 in the following section.

##  V Numerical Simulation

Figure 4: Parameter estimate error for different parameter adaptation gains.

Figure 5: Maximum temperature error for different parameter adaptation gains.

Figure 6: Temperature estimation $L^{2}$ error for different adaptation gains.

To demonstrate the performance of the proposed adaptive observer of (6), (7),
we have run a simulation of an electrosurgical process. The simulation is
performed on a homogeneous domain with size $4\ \text{cm}\times 2\
\text{cm}\times 2\ \text{cm}$ with material properties of porcine loin tissue
from [4], where mass density $\rho=700\ \text{kg}/\text{m}^{3}$, thermal
conductivity $k=0.5934\ \text{W}/\text{kg}\cdot\text{K}$, and isochoric
specific heat $c_{p}=4\ \text{kJ}/\text{kg}\cdot\text{K}$. The heat conduction
is governed by the heat equation in (1), and the entire domain starts at
constant temperature $300\ \text{K}$, with homogeneous Neumann boundary
condition on the top surface and Dirichlet boundary condition on the other
faces. Two cuts are made along the longitudinal axis on the top surface at the
speed of $1\ \text{cm}/\text{s}$. Cut 1 starts from $t=0\ \text{s}$ at
location $(1,0,0)$ and ends at point $(1.75,0,0)$, and cut 2 starts from
$t=1.25\ \text{s}$ at location $(2.25,0,0)$ and ends at point $(3,0,0)$, at
$30\ \text{W}$ output.

The heat source $q$ is a normalized 3D Gaussian distribution defined as
$q(\eta):=q_{0}\exp\left(-\sum_{i=1}^{3}\eta_{i}^{2}/(2\sigma^{2})\right)$,
where $\sigma=1\ \text{mm}$ is the radius of the electrosurgerical probe
needle and $q_{0}$ is a normalization constant that ensures
$\int_{\Omega}q(\eta)\ \text{d}\eta=1$. The domain is discretized spatially
using second-order finite differencing with $\Delta x=0.05\ \text{cm}$, and
temporally using trapezoid integration (Crank–Nicolson) with $\Delta t=0.02\
\text{s}$ and a Gauss–Seidel solver.

First, we run the simulation with the material properties listed above, with
the nominal thermal diffusivity being $a=k/(\rho c_{p})$. To introduce an
initial parameter mismatch, we take $a_{0}=2a$ for the observer system; by
Conjecture 1, we expect the parameter mismatch to be corrected for. As
mentioned in Sec. IV, the ANRA-based thermal diffusivity is only applicable
when the heat source is off, we therefore obtain $\hat{a}(t)$ for
$t\in[0.75,1.25]$ s, i.e., the period when the pen is stationary between both
cuts. In addition, the heat source disturbance $w(t)$ (continuous in time, and
$\gamma$-Hölder in space) is injected into the nominal system of (1) so as to
capture thermographer noise and unmodeled dynamics. The observer gain $G(t)$
is taken to be a constant identity matrix scaled by a factor of $50$. A
video111https://uofi.box.com/s/ghgsv7mmjuf6ogiief1wgslhxroo2lq3 of the
simulated temperature field demonstrates that the adaptation law correctly
detects a change in the thermal diffusivity, with the observer consequently
correcting the estimated temperature field. In addition to this video, Figs.
7(a), 7(b) show evidence of tracking of the decreasing diffusivity, as
indicated by the increased heat spread at $t=1.94$ s.

Fig. 6 shows the error of the adaptive parameter during the simulation for
different values of $L$. We can see that the parameter estimate converges to
the correct value faster as gain $L$ in the adaptive observer increases. In
the case of $L=0.5$, the error goes to zero within one second. However, in the
case of $L=0.1$, there is not enough time for the error to reach zero before
the source is turned back on and parameter adaptation is suspended; given more
time between cuts, this configuration would also have converged.

Referring to the nominal system run, we compute two error metrics
($\infty$-norm and $L^{2}$-norm of the error) over time to demonstrate
convergence of the temperature estimate. In Fig. 6 we can see the effect of
the adaptive observer on the maximum temperature error. We are able to
significantly reduce the error in the second cut even with a small $L$ of 0.1,
and the error goes to zero with $L=0.5$. From Fig. 6, we can see that the
convergence rate is much higher with a larger $L$; towards the end of the
simulation, most of the $L^{2}$ error is due to discretization errors as
evidenced by Fig. 6. In practice, it may not be advantageous to choose
exceedingly large values of $L$, since this may amplify estimation errors due
to sensor noise.

(a) $t=0.54$ s.

(b) $t=1.94$ s.

Figure 7: 3D temperature field in Kelvin at $t=0.54\ \text{s}$ and $t=1.94\
\text{s}$ using adaptive observer (6) with contour slices at planes $y=0$ cm,
and $z=0,-1,-2$ cm.

##  VI Conclusion

In this work, we have considered the observation of a 3D temperature field
arising during an electrosurgical procedure with a moving volumetric heat
source, based only on surface temperature feedback. We have presented a
nonlinear partial differential equation that governs the process, from which
we derived a pointwise noncollocated observer law based on surface
temperature. In this work, we have allowed for a smooth and power-bounded heat
source disturbance to be present across the entire domain, including an
adequate choice of feedback gain. A closed-form lower bound for the feedback
gain has been derived based solely on rudimentary assumptions of the process
noise, as well as the surface temperature feedback. Using this observer law,
asymptotic convergence was shown.

We have also introduced a novel direct parameter estimation technique, which
we termed _attention-based noise-robust averaging_ (ANRA) to enable real-time
parameter adaptation and shown its use. We have demonstrated convergence of
this adaptive observation law through numerical simulation.

## References

  * [1] A. R. Lanfranco, A. E. Castellanos, J. P. Desai, and W. C. Meyers, “Robotic Surgery: A Current Perspective,” _Annals of Surgery_ , vol. 239, no. 1, pp. 14–21, Jan. 2004.
  * [2] D. Palanker, A. Vankov, and P. Jayaraman, “On mechanisms of interaction in electrosurgery,” _New Journal of Physics_ , vol. 10, no. 12, p. 123022, Dec. 2008.
  * [3] H. El-Kebir and J. Bentsman, “PDE-based modeling and non-collocated feedback control of electrosurgical-probe/tissue interaction,” in _2021 American Control Conference_. New Orleans, Louisiana, USA: IEEE, 2021, pp. 4045–4050.
  * [4] H. El-Kebir, Y. Lee, R. Berlin, E. Benedetti, P. C. Giulianotti, L. P. Chamorro, and J. Bentsman, “Online hypermodel-based path planning for feedback control of tissue denaturation in electrosurgical cutting,” in _Proc. of the 11th IFAC Symposium on Biological and Medical Systems_ , B. Benyo, Ed., vol. 54, no. 15. Ghent, Belgium: IFAC, 2021, pp. 448–453.
  * [5] F. Bardati and D. Solimini, “Radiometric sensing of biological layered media,” _Radio Science_ , vol. 18, no. 6, pp. 1393–1401, 1983.
  * [6] A. Erez and A. Shitzer, “Controlled destruction and temperature distributions in biological tissues subjected to monoactive electrocoagulation,” _Journal of Biomechanical Engineering_ , vol. 102, no. 1, pp. 42–49, Feb. 1980.
  * [7] D. Ledwon, A. Sage, J. Juszczyk, M. Rudzki, and P. Badura, “Tomographic reconstruction from planar thermal imaging using convolutional neural network,” _Scientific Reports_ , vol. 12, no. 1, p. 2347, 2022.
  * [8] V. P. Vavilov and X. P. Maldague, “Dynamic thermal tomography: New promise in the IR thermography of solids,” in _Aerospace Sensing_ , J. K. Eklund, Ed., Orlando, FL, USA, 1992, pp. 194–206.
  * [9] J. Toivanen, T. Tarvainen, J. Huttunen, T. Savolainen, H. Orlande, J. Kaipio, and V. Kolehmainen, “3D thermal tomography with experimental measurement data,” _International Journal of Heat and Mass Transfer_ , vol. 78, pp. 1126–1134, 2014.
  * [10] J. Toivanen, T. Tarvainen, J. Huttunen, T. Savolainen, A. Pulkkinen, H. Orlande, J. Kaipio, and V. Kolehmainen, “Thermal tomography utilizing truncated Fourier series approximation of the heat diffusion equation,” _International Journal of Heat and Mass Transfer_ , vol. 108, pp. 860–867, 2017.
  * [11] M. M. Chen, “Computer-aided tomographic thermography,” in _Heat Transfer in Medicine and Biology_ , A. Shitzer and R. C. Eberhart, Eds. Boston, MA, USA: Springer US, 1985, pp. 353–369.
  * [12] S. Singh and R. Melnik, “Thermal ablation of biological tissues in disease treatment: A review of computational models and future directions,” _Electromagnetic Biology and Medicine_ , vol. 39, no. 2, pp. 49–88, Apr. 2020\.
  * [13] P. Holoborodko, “Smooth Noise Robust Differentiators,” http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/, 2008\.
  * [14] ——, “Noise Robust Gradient Operators,” http://www.holoborodko.com/pavel/image-processing/edge-detection/, 2009.