# SRViT: Vision Transformers for Estimating
Radar Reflectivity from Satellite Observations at Scale

Jason Stock  Kyle Hilburn  Imme Ebert-Uphoff  Charles Anderson

###### Abstract

We introduce a transformer-based neural network to generate high-resolution
($3$ km) synthetic radar reflectivity fields at scale from geostationary
satellite imagery. This work aims to enhance short-term convective-scale
forecasts of high-impact weather events and aid in data assimilation for
numerical weather prediction over the United States. Compared to convolutional
approaches, which have limited receptive fields, our results show improved
sharpness and higher accuracy across various composite reflectivity
thresholds. Additional case studies over specific atmospheric phenomena
support our quantitative findings, while a novel attribution method is
introduced to guide domain experts in understanding model outputs.

Machine Learning, Transformers, Deep Learning, Weather and Climate,
Atmospheric Science

##  1 Introduction

Accurate radar is crucial for operational forecasters to monitor and forecast
the progression of high-impact weather events. Given the implications to
public safety and agriculture, among others, this accuracy is key to
protecting life and property. There are two primary uses of radar in weather
forecasting: (1) enabling forecasters to view imagery and decide when to issue
warnings and (2) for imagery to be integrated into numerical weather
prediction (NWP) models to forecast the weather (Gustafsson et al., 2018;
Jones et al., 2020). However, radar is limited to sparsely situated ground
stations, leaving remote areas, mountains, and oceans with poor coverage.
(Figure 7; see Figure 1 in McGovern et al. (2022)). To achieve more accurate
forecasts, it is essential to improve radar coverage both spatially and
temporally, which we aim to address by using observational data.

Figure 1: Slices of an input (col. 1-4) and output prediction with ground
truth (col. 5 and 6, respectively). The real-time observations enable
forecasters to assess/forecast storm patterns at scale. Figure 2: Cropped and
enlarged model output, showing a Northern Plains Derecho in panels (a-c) and
Midwest Squall Lines in panels (d-f). Sample RMSE and $\text{R}^{2}$ values
are shown for each case between the ground truth MRMS (col. 1) and model
output (col. 2 and 3).

The Geostationary Operational Environmental Satellite (GOES) provides
expansive coverage of the contiguous United States (CONUS) with $5$ minute
updates and has been shown to be effective for operational weather forecasting
(Line et al., 2016). With its low-latency and high spatial coverage, we
believe there is a great opportunity in leveraging its geostationary imagers.
In particular, our study focuses on utilizing imagery from the Advanced
Baseline Imager (ABI) (Schmit et al., 2010) and the Geostationary Lightning
Mapper (GLM) (Goodman et al., 2010) onboard GOES-16 to generate high-
resolution synthetic radar as modeled by the Multi-Radar Multi-Sensory (MRMS)
product (Smith et al., 2016) over the CONUS domain (Figure 1).

We use a machine learning approach, a field that has recently received
attention for radar estimation. However, previous work has centered on
incorporating physical indicators from numerical models (Veillette et al.,
2018; Zhu et al., 2023) or modeling localized spatial regions (Hilburn et al.,
2021; Duan et al., 2021; Yu et al., 2023). Despite the influx of high-
resolution weather and climate modeling (Pathak et al., 2022; Mardani et al.,
2024), the transition to high-resolution observational-based modeling is still
in its infancy (Stock et al., 2024). In this work, we hypothesize that
transformer-based models can benefit from contextualizing the synoptic
observations over what is possible with convolutional models. We not only
evaluate our results both quantitatively and through case studies, but we also
introduce an attribution method for transformers to better assist domain
experts.

##  2 Preliminaries

We begin by introducing the dataset in Section 2.1, detailing our data
features and composition, and in Section 2.2 we outline the baseline model
used for making comparisons.

###  2.1 Data Composition

This study uses observational data as input (ABI infrared Channels 7, 9, 13
and GLM) from the Geostationary Operational Environmental Satellite (GOES)-R
Series satellites with target composite radar reflectivity from the Multi-
Radar Multi-Sensor (MRMS) product—complete variable details are provided in
Section A.1. We collect data in part from Hilburn (2023a, b, c) between
2018-2022 that span the contiguous United States (CONUS). As the intent of our
results are to be used with data assimilation, we project the input and target
samples to follow a $3$ km HRRR mass grid, yielding $768\times 1536$-pixel
images ($2304\times 4608$ km) on $6$ hour periods with a $15$ minute refresh
($96$ samples per day).

We restrict data to the warm season (i.e., April-September) where the
strongest radar echos occur. Data are temporally partitioned to the years of
2018-2020 ($47{,}821$) for training, 2021 ($17{,}284$) for validation, and
2022 ($17{,}344$) for testing. Note the number of samples are indicated in
parentheses that altogether total $1.9$ TBs in size. All samples undergo a
preliminary quality control phase to ensure targets have the full range of
radar reflectivity and there are no erroneous samples with large GLM and MRMS
artifacts (additional verification in Section A.2). Additionally, individual
variables are scaled between $[0,1]$ based on their corresponding histograms
to stabilize training. When unstandardized, MRMS and model output values scale
between $[0,60]$ dBZ (unit of reflectivity).

###  2.2 Baseline Model Comparison

The GOES Radar Estimation via Machine Learning to Inform NWP (GREMLIN) model
(Hilburn et al., 2021) uses a fully-convolutional encoder/decoder architecture
(without skip connections) that we denote as “UNet” in this study. Trained
previously on smaller $256\times 256$-pixel images and regionally filtered by
storm reports (tornado, hail, and wind events), this architecture serves as
the baseline for satellite to radar emulation. Using the same hyperparameters
as in the original work, we retrain the UNet and both evaluate and compare its
performance to our model using data from the entire CONUS domain.

Figure 3: Categorical metrics at varying composite reflectivity thresholds for
SRViT and the baseline UNet.

Figure 4: Kernel density estimation (KDE) of the mean gradient magnitude of
composite reflectivity over all test samples. The dashed line represents the
mean with standard deviation.

##  3 Methodology

Our proposed model extends the vision transformer (Vaswani et al., 2017;
Dosovitskiy et al., 2020) and is designed specifically for image-to-image
translation. This is motivated by wanting to address the limitations of
restricted receptive fields and empirical evidence for blurriness in tandem
with pixel-wise loss functions (Ravuri et al., 2021; Ayzel et al., 2020) with
purely convolutional-based methods. By using self-attention, we are able to
contextualize the long-range dependencies and synoptic information of global
weather patterns that influence local precipitation events.

The design specifics of our Satellite to Radar Vision Transformer (SRViT), as
detailed in Appendix B, is inspired by prior global weather forecasting
models. In particular, work from Pathak et al. (2022); Bi et al. (2022); Chen
et al. (2023) that all use transformer-based architectures and are shown to
outperform graph- and convolutional-based approaches while approaching the
performance of NWP (Lam et al., 2022; Rasp & Thuerey, 2021). However, in this
work we simplify the architecture and training procedure.

We use the weighted loss from Hilburn et al. (2021) to balance the rare but
high radar reflectivity values with the small but common values. Specifically,
denoted as

$\displaystyle\mathcal{L}_{e}=\frac{1}{m}\sum_{i=1}^{m}{\exp{(w_{0}t_{i}^{w_{1}})}\cdot(y_{i}-t_{i})^{2}},$
(1)

where $t$ and $y$ are the ground truth and predicted values, $m$ is the number
of pixels, and $w_{0}=5$ and $w_{1}=4$ are found by optimizing the categorical
bias of model performance.

We train our $0.65$M parameter model on a cluster of $8\times$ $40$ GB NVIDIA
A100 GPUs for a maximum of $300$ epochs with AdamW. End-to-end training takes
about $100$ hours wall-clock time, and due to high compute and limited
accessibility, we perform minimal hyperparameter tuning (see Table 2 for final
values) and report results of a single trial.

##  4 Experimental Results

Primary results and comparisons with the baseline model are discussed in
Section 4.1, and in Section 4.2 we investigate individual case studies of
severe weather events (e.g., Figure 2). Additional examples of model output in
Appendix C.

###  4.1 Main Findings

We evaluate our results with three unique categories of methods, including:
(a) standard metrics: root-mean-squared error (RMSE) and coefficient of
determination ($\text{R}^{2}$); (b) categorical metrics: probability of
detection (POD), false alarm ratio (FAR), and critical success index (CSI) at
different composite reflectivity thresholds; and (c) sharpness quantification
with the mean magnitude of image gradients.

#### Model Performance

SRViT overall has better performance on the standard metrics, achieving an
$\text{RMSE}=3.09\text{ dBZ and r}^{2}=0.572$, compared to the UNet’s
$\text{RMSE}=3.21\text{ dBZ and r}^{2}=0.488$ (Table 2). These statistics
capture general pixel-wise improvement, but are biased toward low coverage and
zero-valued pixels. To improve our intuition of model performance, we evaluate
the results on a linear interval where composite reflectivity is greater than
predefined threshold levels. This more accurately displays where the
improvements are captured.

Notably, in Figure 4, we find SRViT exhibits the best performance at low- to
mid-value thresholds. Relative to the UNet, the greatest improvements
represented by FAR, CSI, and RMSE are between $[5,40)$ dBZ. Specifically,
there is a reduction of $0.96$ dBZ in RMSE for reflectivity greater than $35$
dBZ, while FAR increases by $0.08$ at a threshold of $25$ dBZ. However, the
POD of each network is consistent across thresholds until a divergence occurs
at $40$ dBZ. At these higher thresholds ($>40$ dBZ), the UNet displays
marginally better performance across all the categorical indicators, which may
be attributed to specific test cases.

#### Sharpness Evaluation

Qualitative sharpness can be defined by distinct, and sharp edge boundaries of
neighboring pixels. An observation we find the lack of with the boundaries of
composite reflectivity from the UNet. To quantify this observation, we compute
the mean of the gradient magnitude for each output and compare the sample
distribution among models. Specifically, we take the spatial mean after
convolving a Sobel filter, i.e,
$g=\frac{1}{m}\sum_{i=1}^{m}(G_{x_{i}}^{2}+G_{y_{i}}^{2})^{\frac{1}{2}}$,
where $x$ and $y$ are image directions, for each sample. Figure 4 shows the
kernel density estimation (KDE) plots of $g$ across all test samples for each
model.

Both networks exhibit positively skewed distributions with lower mean values
compared to MRMS. A direct comparison using Welch’s independent samples
t-test, show that SRViT demonstrates a statistically significant improvement
in sharpness over the UNet ($\text{t-statistic}=27.71\text{, p-value}<0.001$).
These finding coincide with the qualitative assessments of observations (e.g.,
Figure 2) and suggest that our method enhances the sharpness of composite
reflectivity. However, both approaches still fall short of fully capturing the
sharpness within the ground truth MRMS data.

###  4.2 Case Studies

By enlarging model output, we study the results of our approach over a small
spatial region for two distinct severe weather events and compare the
predictions with MRMS and the UNet. In Figure 2, the top row shows a Northern
Plains Derecho and the bottom row shows Midwest Squall Lines. In both cases,
SRViT has an improvement to the standard metrics (lower RMSE, higher
$\text{R}^{2}$) relative to the UNet. A qualitative assessment has composite
reflectivity boundaries with precise transitions, validating our findings of
generating sharper output. The UNet also appears to overestimate for higher
thresholds, whereas SRViT more accurately captures low- to mid-level
thresholds but underestimates composite reflectivity $>50$ dBZ, e.g., within
the lower squall line over Southern Ohio in the bottom row of Figure 2. These
observations are consistent with the quantitative metrics we report in Section
4.1.

##  5 Interpreting Model Attribution

Some prior work aims to interpret the attention matrix of transformer models
as importance or relevancy scores (Kovaleva et al., 2019; Clark et al., 2019;
Carion et al., 2020). Indeed, these weights capture the interactions between
all tokens, but faithful interpretations are challenged by not considering (a)
the token vector magnitudes and high dimensionality and (b) the interplay
between tokens across intermediate hidden layers (Brunner et al., 2019;
Wiegreffe & Pinter, 2019; Serrano & Smith, 2019).

We address these issues by introducing Token (Re)Distribution: an attribution
method that captures the interaction of all other tokens, as learnt by the
network, for a particular token. This approach allows for interpretations of
how the value of input tokens are redistributed amongst the others. Based on
prior work (Hao et al., 2021; Dhamdhere et al., 2018; Chefer et al., 2021)
that use the gradient of intermediate layers and self-attention features, we
summarize gradient information across token embeddings to highlight the
relevance of individual tokens (see details in Appendix E).

An example attribution map with model output of Tropic Storm Colin is shown in
Figure 5. Each panel shows the result of a different token, from the same
sample, with the normalized gradient magnitude. In Figure 5a the token of
interest, centered on higher radiance values, shows surrounding tokens within
the clouds contributing to its value. Similarly, in Figure 5b, the token on
the storm edge uses the values of others along this edge. Both cases are
meteorological supported as it relates to the storm structure and extent of
precipitation. Altogether, this method is a step toward guiding domain
experts’ interpretations of model output.

Figure 5: Sample min-max normalized gradient magnitude of Token
(Re)Distribution for the (red; separate in panels (a) and (b)) token of
interest overlaid on GOES-16 ABI Channel 7.

##  6 Conclusion

Our analyses demonstrate that SRViT improves low- and mid-value estimates of
composite reflectivity and produces sharper predictions compared to
convolutional approaches. Case studies validate the model’s capability to
accurately capture complex atmospheric phenomena both locally and across
larger domains. Additionally, our gradient-based attribution method provides
valuable meteorological insights and is applicable to various transformer
architectures.

Despite these advancements, some limitations remain. The UNet outperforms
SRViT at high-value estimates of composite reflectivity, and neither model
fully achieves the sharpness of MRMS. Future work could address these gaps by
exploring probabilistic diffusion models to better capture uncertainty and
sharpness while also developing methods to estimate 3D radar reflectivity
fields to better inform NWP.

Code available on GitHub: github.com/stockeh/srvit

## Acknowledgments

This work is supported by NSF Grant No. 2019758, AI Institute for Research on
Trustworthy AI in Weather, Climate, and Coastal Oceanography (AI2ES) as well
as the GOES-R Program under grant NA19OAR4320073. A special thanks to the
University at Albany for allowing us to use their xCITE A100 GPU compute
cluster for experiments and model training.

## References

  * Ayzel et al. (2020) Ayzel, G., Scheffer, T., and Heistermann, M.  Rainnet v1. 0: a convolutional neural network for radar-based precipitation nowcasting.  _Geoscientific Model Development_ , 13(6):2631–2644, 2020.
  * Bi et al. (2022) Bi, K., Xie, L., Zhang, H., Chen, X., Gu, X., and Tian, Q.  Pangu-weather: A 3d high-resolution model for fast and accurate global weather forecast.  _arXiv preprint arXiv:2211.02556_ , 2022.
  * Brunner et al. (2019) Brunner, G., Liu, Y., Pascual, D., Richter, O., Ciaramita, M., and Wattenhofer, R.  On identifiability in transformers.  _arXiv preprint arXiv:1908.04211_ , 2019.
  * Carion et al. (2020) Carion, N., Massa, F., Synnaeve, G., Usunier, N., Kirillov, A., and Zagoruyko, S.  End-to-end object detection with transformers.  In _European conference on computer vision_ , pp. 213–229. Springer, 2020.
  * Chefer et al. (2021) Chefer, H., Gur, S., and Wolf, L.  Transformer interpretability beyond attention visualization.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 782–791, 2021.
  * Chen et al. (2023) Chen, K., Han, T., Gong, J., Bai, L., Ling, F., Luo, J.-J., Chen, X., Ma, L., Zhang, T., Su, R., et al.  Fengwu: Pushing the skillful global medium-range weather forecast beyond 10 days lead.  _arXiv preprint arXiv:2304.02948_ , 2023.
  * Clark et al. (2019) Clark, K., Khandelwal, U., Levy, O., and Manning, C. D.  What does bert look at? an analysis of bert’s attention.  _arXiv preprint arXiv:1906.04341_ , 2019.
  * Dhamdhere et al. (2018) Dhamdhere, K., Sundararajan, M., and Yan, Q.  How important is a neuron?  _arXiv preprint arXiv:1805.12233_ , 2018.
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Duan et al. (2021) Duan, M., Xia, J., Yan, Z., Han, L., Zhang, L., Xia, H., and Yu, S.  Reconstruction of the radar reflectivity of convective storms based on deep learning and himawari-8 observations.  _Remote Sensing_ , 13(16):3330, 2021.
  * Goodman et al. (2010) Goodman, S. J., Mach, D., Koshak, W., and Blakeslee, R.  Glm lightning cluster-filter algorithm (lcfa) algorithm theoretical basis document (atbd).  _NOAA NESDIS STAR Document_ , pp. 77, 2010.
  * Goodman et al. (2013) Goodman, S. J., Blakeslee, R. J., Koshak, W. J., Mach, D., Bailey, J., Buechler, D., Carey, L., Schultz, C., Bateman, M., McCaul Jr, E., et al.  The goes-r geostationary lightning mapper (glm).  _Atmospheric research_ , 125:34–49, 2013.
  * Gustafsson et al. (2018) Gustafsson, N., Janjić, T., Schraff, C., Leuenberger, D., Weissmann, M., Reich, H., Brousseau, P., Montmerle, T., Wattrelot, E., Bučánek, A., et al.  Survey of data assimilation methods for convective-scale numerical weather prediction at operational centres.  _Quarterly Journal of the Royal Meteorological Society_ , 144(713):1218–1256, 2018.
  * Hao et al. (2021) Hao, Y., Dong, L., Wei, F., and Xu, K.  Self-attention attribution: Interpreting information interactions inside transformer.  In _The Thirty-Fifth AAAI Conference on Artificial Intelligence_. AAAI Press, 2021.  URL https://arxiv.org/pdf/2004.11207.pdf.
  * Hilburn (2023a) Hilburn, K.  Gremlin conus3 dataset for 2020, 2023a.  Dryad, Dataset, https://doi.org/10.5061/dryad.h9w0vt4nq.
  * Hilburn (2023b) Hilburn, K.  Gremlin conus3 dataset for 2021, 2023b.  Dryad, Dataset, https://doi.org/10.5061/dryad.zs7h44jf2.
  * Hilburn (2023c) Hilburn, K.  Gremlin conus3 dataset for 2022, 2023c.  Dryad, Dataset, https://doi.org/10.5061/dryad.2jm63xstt.
  * Hilburn et al. (2021) Hilburn, K. A., Ebert-Uphoff, I., and Miller, S. D.  Development and interpretation of a neural-network-based synthetic radar reflectivity estimator using goes-r satellite observations.  _Journal of Applied Meteorology and Climatology_ , 60(1):3 – 21, 2021.  doi: https://doi.org/10.1175/JAMC-D-20-0084.1.  URL https://journals.ametsoc.org/view/journals/apme/60/1/jamc-d-20-0084.1.xml.
  * Jones et al. (2020) Jones, T. A., Skinner, P., Yussouf, N., Knopfmeier, K., Reinhart, A., Wang, X., Bedka, K., Smith Jr, W., and Palikonda, R.  Assimilation of goes-16 radiances and retrievals into the warn-on-forecast system.  _Monthly Weather Review_ , 148(5):1829–1859, 2020.
  * Kovaleva et al. (2019) Kovaleva, O., Romanov, A., Rogers, A., and Rumshisky, A.  Revealing the dark secrets of bert.  _arXiv preprint arXiv:1908.08593_ , 2019.
  * Lam et al. (2022) Lam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Pritzel, A., Ravuri, S., Ewalds, T., Alet, F., Eaton-Rosen, Z., et al.  Graphcast: Learning skillful medium-range global weather forecasting.  _arXiv preprint arXiv:2212.12794_ , 2022.
  * Latto & Cangialosi (2022) Latto, A. S. and Cangialosi, J. P.  National hurricane center tropical cyclone report: Tropical storm colin.  _Miami: National Hurricane Center_ , 2022.
  * Line et al. (2016) Line, W. E., Schmit, T. J., Lindsey, D. T., and Goodman, S. J.  Use of geostationary super rapid scan satellite imagery by the storm prediction center.  _Weather and Forecasting_ , 31(2):483–494, 2016.
  * Mardani et al. (2024) Mardani, M., Brenowitz, N., Cohen, Y., Pathak, J., Chen, C.-Y., Liu, C.-C., Vahdat, A., Kashinath, K., Kautz, J., and Pritchard, M.  Residual diffusion modeling for km-scale atmospheric downscaling.  _arXiv preprint arXiv:2309.15214_ , 2024.
  * McGovern et al. (2022) McGovern, A., Ebert-Uphoff, I., Gagne, D. J., and Bostrom, A.  Why we need to focus on developing ethical, responsible, and trustworthy artificial intelligence approaches for environmental science.  _Environmental Data Science_ , 1:e6, 2022.
  * Pathak et al. (2022) Pathak, J., Subramanian, S., Harrington, P., Raja, S., Chattopadhyay, A., Mardani, M., Kurth, T., Hall, D., Li, Z., Azizzadenesheli, K., et al.  Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators.  _arXiv preprint arXiv:2202.11214_ , 2022.
  * Rasp & Thuerey (2021) Rasp, S. and Thuerey, N.  Data-driven medium-range weather prediction with a resnet pretrained on climate simulations: A new model for weatherbench.  _Journal of Advances in Modeling Earth Systems_ , 13(2):e2020MS002405, 2021.
  * Ravuri et al. (2021) Ravuri, S., Lenc, K., Willson, M., Kangin, D., Lam, R., Mirowski, P., Fitzsimons, M., Athanassiadou, M., Kashem, S., Madge, S., et al.  Skilful precipitation nowcasting using deep generative models of radar.  _Nature_ , 597(7878):672–677, 2021.
  * Schmit et al. (2010) Schmit, T., Gunshor, M., Fu, G., Rink, T., Bah, K., and Wolf, W.  Goes-r advanced baseline imager (abi) algorithm theoretical basis document for: Cloud and moisture imagery product (cmip), version 2.3.  _NOAA NESDIS STAR Document_ , pp. 66, 2010.
  * Serrano & Smith (2019) Serrano, S. and Smith, N. A.  Is attention interpretable?  _arXiv preprint arXiv:1906.03731_ , 2019.
  * Smith et al. (2016) Smith, T. M., Lakshmanan, V., Stumpf, G. J., Ortega, K. L., Hondl, K., Cooper, K., Calhoun, K. M., Kingfield, D. M., Manross, K. L., Toomey, R., et al.  Multi-radar multi-sensor (mrms) severe weather and aviation products: Initial operating capabilities.  _Bulletin of the American Meteorological Society_ , 97(9):1617–1630, 2016.
  * Stock et al. (2024) Stock, J., Pathak, J., Cohen, Y., Pritchard, M., Garg, P., Durran, D., Mardani, M., and Brenowitz, N.  Diffobs: Generative diffusion for global forecasting of satellite observations.  _arXiv preprint arXiv:2404.06517_ , 2024.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Veillette et al. (2018) Veillette, M. S., Hassey, E. P., Mattioli, C. J., Iskenderian, H., and Lamey, P. M.  Creating synthetic radar imagery using convolutional neural networks.  _Journal of Atmospheric and Oceanic Technology_ , 35(12):2323–2338, 2018.
  * Wiegreffe & Pinter (2019) Wiegreffe, S. and Pinter, Y.  Attention is not not explanation.  _arXiv preprint arXiv:1908.04626_ , 2019.
  * Yu et al. (2023) Yu, X., Lou, X., Yan, Y., Yan, Z., Cheng, W., Wang, Z., Zhao, D., and Xia, J.  Radar echo reconstruction in oceanic area via deep learning of satellite data.  _Remote Sensing_ , 15(12):3065, 2023.
  * Zhu et al. (2023) Zhu, M., Liao, Q., Wu, L., Zhang, S., Wang, Z., Pan, X., Wu, Q., Wang, Y., and Su, D.  Multiscale representation of radar echo data retrieved through deep learning from numerical model simulations and satellite images.  _Remote Sensing_ , 15(14):3466, 2023.

##  Appendix A Dataset Details

Section A.1 overviews the critical data variables and properties unique to
this study; however, we refer the reader to Hilburn et al. (2021) for more
complete data specifications. We also provide insights to the spatial coverage
in Section A.2.

###  A.1 Variable Description

#### Advanced Baseline Imager (ABI)

We use Level-L1b radiances from GOES-16 ABI (Schmit et al., 2010). The imagery
has a nominal spatial resolution of $2$ km, enabling the capture of fine-
grained details, while ensuring a temporal refresh rate of $5$ minutes for
timely data updates. To account for both day-night scenarios, we concentrate
on the use of the infrared bands, specifically: Channel 7
($3.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$; Shortwave Window), Channel 9
($6.9\text{\,}\mathrm{\SIUnitSymbolMicro m}$; Upper-Level Water Vapor), and
Channel 13 ($10.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$; Clean IR Longwave
Window).

#### Geostationary Lightning Mapper (GLM)

Accompanying the satellite imagery are real-time lightning observations from
the GOES-R GLM (Goodman et al., 2010, 2013). This instrument is a single-
channel, near-infrared transient detector that monitors total lightning with a
uniform spatial resolution of approximately $10$ km. Lightning itself is
particularly useful for generating synthetic radar fields as their locations
are often associated with strong updrafts within convective environments.

#### Multi-Radar Multi-Sensor (MRMS)

We utilize quality-controlled composite reflectivity from the MRMS product
(Smith et al., 2016) as our target dataset. This data combines information
from different radar networks, surface observations, numerical weather
prediction (NWP) models, and climatology. Through a sophisticated integration
process, high-resolution mosaics are generated, offering precise
spatiotemporal resolution. However, it is important to note that the coverage
of the product is both limited in the vertical and spatial domain due to beam
blockage and radar placement.

###  A.2 Data Coverage and Distribution

Figure 7 illustrates the radar coverage from the training data, presenting a
spatial map of cumulative echos. To represent the spatial distribution, we
binarize the data by assigning a value of $1$ to composite reflectivity values
greater than $0$, which are then summed across all samples. We then clip the
result based on the $99{\text{-th}}$ percentile to eliminate visual artifacts.
Notably, the coverage is lowest over the western United States, with a bias
towards physical radar stations. The western region, affected by orographic
enhancement, is further impeded by mountains (i.e., subject to beam blockage),
resulting in predominantly high-level atmospheric echoes that largely miss the
stronger echoes near the surface. In the most extreme cases, precipitation may
go completely undetected by radars.

In order to disentangle the sample specific coverage from the cumulative
spatial representation, we compute the sample frequency of nonzero pixels with
its cumulative density function (CDF). Figure 7 displays the result on the
training data, showing a histogram with a slight positive skew and mean
centered around $9.2\%$, indicating moderate coverage within samples relative
to the CONUS domain. Moreover, a substantial proportion of the training
samples have nonzero pixel percentages above the median value. Importantly,
this indicates that all samples have some level of coverage, irrespective of
their spatial extent, without exceeding into high, erroneous values.

Figure 6: Spatial coverage of MRMS from the training data.

Figure 7: Sample histogram and cumulative density function of nonzero pixels
present in the training data.

Hyperparameter | Value
---|---
Global Batch Size | 16
Learning Rate $\eta$ | $1e-4$
Epochs (patience) | 300 (50)
Patch Size $(p\times p)$ | $(12\times 12)$
Model Depth $L$ | 6
Heads (per block) | 12
Model Dimension $d$ | 256
Inner Dimension $v$ | 64
FCN Dimension $m$ | 512
Conv Hiddens (filters) | $[32,16]$
Table 1: SRViT hyperparameters ($643{,}105$ parameters).

Model |  $\downarrow$ RMSE (dBZ) | $\uparrow\text{R}^{2}$ |  $\uparrow$ Sharpness ($g$)
---|---|---|---
MRMS | – | – |  0.48 $\pm$ 0.16
UNet | 3.21 | 0.488 |  0.21 $\pm$ 0.09
Base-ViT | 3.05 | 0.487 |  0.21 $\pm$ 0.09
SRViT | 3.09 | 0.572 |  0.24 $\pm$ 0.11
Table 2: Standard and sharpness metrics across models.

##  Appendix B Architectural Details

Consider an image $I\in\mathbb{R}^{c\times h\times w}$ that is partitioned
into equally divisible patches (synonymous with tokens) of size $p$ as
$x=\\{x_{0},\dots,x_{n}\\}$, where $x_{i}\in\mathbb{R}^{d_{in}}$ and
$n=(w/p)\cdot(h/p)$ with $d_{in}=p^{2}c$. Once partitioned, the set of tokens
$x$ are linearly projected to dimension $d<d_{in}$ and we add a standard sine-
cosine positional encoding. Let matrix $\mathbf{X}^{l}\in\mathbb{R}^{n\times
d}$ be the new row-wise concatenation of the tokens. A typical transformer
block $\phi$ at layer $l$ processes the set of tokens with multi-head self
attention (MSA) and a point-wise fully-connected network (FCN) as,

$\displaystyle\phi(\mathbf{X})$
$\displaystyle=\text{FCN}(\text{MSA}(\mathbf{X}))\quad\text{such that}$ (2)
$\displaystyle\text{MSA}(\mathbf{X})$
$\displaystyle=[\mathbf{O}_{1},\mathbf{O}_{2},\dots,\mathbf{O}_{h}]\mathbf{W^{O}},$
(3)

where $h$ is the number of heads, $\mathbf{W^{O}}\in\mathbb{R}^{hv\times d}$
are trainable weights, $[\cdot]$ is the column-wise concatenation, and
$\mathbf{O}_{i}\in\mathbb{R}^{n\times v}$ is the output of the $i$-th
attention head with latent dimension $v<d$. We compute each head as,

$\displaystyle\mathbf{O}_{i}$
$\displaystyle=\mathbf{A}_{i}\mathbf{V}_{i}\quad\text{such that}$ (4)
$\displaystyle\mathbf{A}_{i}$
$\displaystyle=\text{softmax}\big{(}\mathbf{Q}_{i}\mathbf{K}_{i}^{\top}/\sqrt{d}\big{)}\in\mathbb{R}^{n\times
n}.$ (5)

The queries, $\mathbf{Q}_{i}$, keys, $\mathbf{K}_{i}$, and values,
$\mathbf{V}_{i}$ are found via a linear projection of $\mathbf{X}$ by,

$\mathbf{Q}_{i}=\mathbf{X}\mathbf{W}^{Q}_{i},\quad\mathbf{K}_{i}=\mathbf{X}\mathbf{W}^{K}_{i},\quad\mathbf{V}_{i}=\mathbf{X}\mathbf{W}^{V}_{i},$
(6)

with trainable weight matrices
$\mathbf{W}^{Q}_{i},\mathbf{W}^{K}_{i},\mathbf{W}^{V}_{i}\in\mathbb{R}^{d\times
v}$. The FCN block at layer $l$ is a two layer network separated by the
Gaussian Error Linear Unit (GELU) activation, $\delta$, and dropout (with a
default $p=0.2$) that takes as input the layer normalized output of the MSA,
$\bar{\mathbf{X}}=\text{LN}(\mathbf{X})$. This block is computed as,

$\text{FCN}(\bar{\mathbf{X}})=\delta(\bar{\mathbf{X}}\mathbf{W}^{R})\mathbf{W}^{S},$
(7)

where $\mathbf{W}^{R}\in\mathbb{R}^{d\times m}$ and
$\mathbf{W}^{S}\in\mathbb{R}^{m\times d}$ such that $m>d$.

Each transformer block has independent weight matrices and yields a new set of
tokens of the same dimension given by
$\phi:\mathbf{X}^{l}\rightarrow\mathbf{X}^{l+1}$. After $L$ transformations,
we introduce a linear decoding block and reshaping operation $f$ to
reconstruct the intermediate output. This is denoted by,

$z_{0}=f(\mathbf{X}^{L})=\text{reshape}(\mathbf{X}^{L}\mathbf{W}^{F},(c\times
h\times w)),$ (8)

where $\mathbf{W}^{F}\in\mathbb{R}^{d\times d_{in}}$ and $z_{0}$ is the same
dimension as the input image $I$. The output of this layer is a close
approximation to the output, but as evidence in Appendix D, we show it to be
inaccurate. Therefore, we include convolutional layers following $f$ to
effectively smooth the boundaries of the decoded tokens via weight sharing. We
represent the $N$ subsequent convolutional layers as,

$z_{i}=\text{ReLU}(\mathbf{W}^{C}_{i}*z_{i-1}),\text{ for }i=1\text{ to }N,$
(9)

where $\ast$ denotes the convolution operation, $\mathbf{W}^{C}_{i}$ are the
trainable weights at layer $i$, and ReLU is applied element-wise. The model’s
final output, $y$, when $i=N$ has the same spatial dimension as $I$ with a
single output channel.

##  Appendix C Additional Model Output

Figure 8: Example model output, showing the temporal progression of Tropical
Storm Colin as it compares MRMS (top row) and SRViT (bottom row). Sample RMSE
and $\text{R}^{2}$ values are shown for each timestep (computed between
individual panel columns).

Figure 8 illustrates the progression of Tropical Storm Colin on 12 hour
intervals from the output of SRViT on an Equidistant Cylindrical map
projection. The short-lived storm had a peak intensity of $35$ knots from 2330
UTC July 1, 2023 to 1200 UTC July 2, 2023 over the western Atlantic.
Accompanied with deep convection, Colin brought at most $19.28$ cm of rainfall
to Wadmalaw Island, South Carolina (Latto & Cangialosi, 2022). SRViT
adequately captures Colin’s structure (Figure 8e) while accurately
representing inland scattered convection (Figure 8g,h). This observation
demonstrates its ability to handle complex weather patterns over the entire
CONUS domain.

##  Appendix D Ablation Study

To better understand our network design, we ablate the key components and
discuss the qualitative and quantitative implications. Following the
transformer blocks in SRViT are convolutional layers that effectively smooth
the boundaries of decoded tokens. By training the base transformer model with
the same hyperparameters and training procedure but without convolutions, we
can study this observation.

We find quantitative statistics of the base model having an
$\text{RMSE}=3.05\text{ dBZ and r}^{2}=0.487$ (see Table 2). This shows a
lower pixel-wise error, which means its predictions are, on average, closer to
MRMS. However, the lower $\text{R}^{2}$ suggests that the use of convolutions
better explains a larger proportion of the variance and captures the
underlying patterns more effectively. When assessing sharpness, it becomes
evident that the base model, having a mean gradient magnitude of $0.21\pm
0.09$, fails to capture the subtle transitions of composite reflectivity
boundaries.

Figure 9: Observational “patchiness” most evident within the yellow box across
model architectures.

Through qualitative observations (Figure 9), it is evident that the base model
produces output with noticeable tiling or “patchiness”. We suspect this
artifact arises from the process of decoding each token individually, which
negatively impacts the natural transition and coherence of the output. By
contrast, incorporating additional convolutional layers proves to be effective
in capturing intricate details and avoiding the tiling issues, making it a
better choice for generating synthetic radar.

##  Appendix E Details of Token (Re)Distribution

The activation of scaled dot-product attention computes the matrix
$\mathbf{A}\in\mathbb{R}^{n\times n}$ that represents the reweighing of our
values from the combined queries and keys (see Equations 4 and 5). Recall that
each transformer block yields an activation matrix as given by
$\phi:\mathbf{X}^{l}\rightarrow\mathbf{X}^{l+1}\in\mathbb{R}^{n\times d}$. Let
the input token embedding prior to adding the positional encoding be
$\mathbf{X}$ and the output of the first block be $\mathbf{Z}$ (this can be
for any block following the rules for backpropagation). The output matrix
$\mathbf{U}\in\mathbb{R}^{n\times n}$ has for the $i$-th intermediate token
(as a row-vector, $\mathbf{z}$), the sensitivity to all input tokens in the
$j$-th column. By differentiating $\partial\mathbf{z}/\partial\mathbf{X}$, we
encode the partial derivatives in a third-order tensor,
$\mathbf{D}^{\mathbf{z}}\in\mathbb{R}^{d\times n\times d}$, indexed by
$\mathbf{D}^{\mathbf{z}}_{i,j,k}=\partial\mathbf{z}_{i}/\partial\mathbf{X}_{jk}$.
In practice, we compute this efficiently as the vector-Jacobian product with
the row-vector $\mathbf{1}=[1_{1},\dots,1_{d}]$ that reduces the tensor to
$\mathbf{D}^{\mathbf{z}}\in\mathbb{R}^{n\times d}$. We summarize the
computation with all $n$ row-vectors in $\mathbf{Z}$ as,

$\displaystyle\mathbf{U}$
$\displaystyle=[f(\mathbf{Z}_{1})_{1},\dots,f(\mathbf{Z}_{n})_{n}]\quad\text{such
that}$ (10) $\displaystyle f(\mathbf{z})_{i}$
$\displaystyle=\Big{\lvert}\sum\nolimits_{j=1}^{d}\mathbf{D}^{\mathbf{z}}_{ij}\Big{\rvert}\in\mathbb{R}^{n}.$
(11)

The function $f$ reduces the $d$-dimensional gradients of each token’s
sensitivity for all $n$ input tokens. This provides a holistic view for how
the change in the input token, over its entire embedding, affects the
intermediate token; a large magnitude is highly responsive to small changes.
Thus, we define a redistribution of those input tokens, as a result of self-
attention, to the value of an intermediate token. With multiple transformer
blocks, we can take the mean over all blocks for a more complete view.
Thereafter, visualization of this mean are made using data indexed at a given
row in $\mathbf{U}$.


# An elastically stabilized spherical invagination

Xiaoyu Zheng1, Tianyi Guo2 , Peter Palffy-Muhoray1,2
_1 Department of Mathematical Sciences, Kent State University, OH, USA _
_2 Advanced Materials and Liquid Crystal Institute, Kent State University, OH,
USA_ Email<EMAIL_ADDRESS><EMAIL_ADDRESS>author Email:
<EMAIL_ADDRESS>
###### Abstract

Invaginations are partial enclosures formed by surfaces. Typically formed by
biological membranes; they abound in nature. In this paper, we consider
fundamentally different structures: elastically stabilized invaginations.
Focusing on spherical invaginations formed by elastic membranes, we carried
out experiments and mathematical modeling to understand the stress and strain
fields underlying stable structures. Friction plays a key role in
stabilization, and consequently the required force balance is an inequality.
Using a novel scheme, we were able to find stable solutions of the balance
equations for different models of elasticity, with reasonable agreement with
experiments.

Keywords: hyperelasticity, invagination, friction

##  1 Introduction

Invaginations are cavities with narrow openings; they are partial enclosures
formed by surfaces. They abound in nature. In the process of pinocytosis [1],
a cell absorbs particles - nutrients - by surrounding them with its membrane,
first forming an invagination partially enclosing the particle and
subsequently completely enclosing and engulfing it. Such invaginations are
typically created by membranes and not by the passive partially enclosed
object. Here we consider a different process, where the deformation of an
elastic membrane is the result of forces exerted by a solid body resulting in
an invagination entrapping the deforming body. An elastically stabilized
deformation of an elastic membrane surrounding a rigid disk has already been
considered [2]. Motivation for this work is to explore the effects of
geometry. Here we study a spherical invagination in an elastic membrane
surrounding a rigid sphere, as shown in Fig. 1, both experimentally and
theoretically. We first describe our experiments in producing and
characterizing spherical invaginations. We next discuss simple mathematical
models and numerical simulations. We close by summarizing our findings.

Figure 1: An elastically stabilized spherical invagination formed by a dental
dam and a table tennis ball.

##  2 Experiments

###  2.1 The membrane and the sphere

The elastic membranes used in our experiments were thin natural latex sheets
sold as powder-free Crosstex 19202 dental dams. Their dimensions are
$5^{\prime\prime}\times 5^{\prime\prime}\times 150\mu m$. Their elastic
modulus is $2MPa$111 The modulus and thickness of the membranes were measured
by us.. The maximum stretch we were able to achieve in one direction without
damaging the membrane is $6$. The spheres used in the experiments discussed
were regulation size $40mm$ dia. table tennis balls. Since friction between
the ball and the membrane plays a key role in stabilizing the invagination, we
attempted to measure both the static $\mu_{s}$ and the kinetic $\mu_{k}$
coefficients of friction between the membrane and the ball. For the
measurements, we used narrow strips ($15cm\times 5mm$) of the membrane. For
the static measurements, we attached a mass at each end of the strip supported
by the ball. We increased the mass at one end until the membrane started slip.
The friction coefficient $\mu_{s}$ can be obtained from the ratio of the
masses via the Euler capstan equation [3],

$T_{load}=T_{hold}\exp(\mu\theta),$ (1)

where $T$ is the tension, $\theta$ is the subtended angle. For the kinetic
measurements, we used a lathe to rotate the table tennis ball while it
supported the strip. One end of the strip was fixed and held horizontal while
a mass was suspended from the other. We measured the strain of the strip near
the fixed end while ball was rotating, $\mu_{k}$ can be obtained from the
ratio of the force producing the strain and the weight of the mass again via
Eq. (1). Our measurements gave $\mu_{s}=0.99$ and $\mu_{k}=0.82$. We were
unable to find results for the kinetic coefficient in the literature, but our
value for the static coefficient is in the range $0.86<\mu_{k}<2.18$ reported
in [4].

###  2.2 Forming the invagination

The invagination by the membrane was realized using two distinct methods.

1\. Our original method consisted of stretching the central portion of the
membrane in plane, by hand, essentially producing an extremely highly
stretched horizontal circular region, approximately $5^{\prime\prime}\times
5^{\prime\prime}$, in the center of the original sheet. This was then placed
above the table tennis ball, which was supported by the neck of a glass bottle
with $1^{\prime\prime}$ OD opening. The stretched membrane was lowered onto
the ball, maintaining cylindrical symmetry and maintaining the radial stretch
as the edges were further lowered below and beneath the ball; necking of the
stretched membrane below the ball could be clearly observed. The tension was
then slowly released, with the membrane contracting and slipping on the ball,
coming to rest and forming the invagination shown in Fig. 1.

2\. Here the unstretched membrane was placed on one end of a
$4^{\prime\prime}$ high $4^{\prime\prime}$ OD PVC tube. The membrane was held
in place by a $1/2^{\prime\prime}$ wide rubber band pressing the edges of the
membrane to the tube. With the tube vertical and the membrane at the bottom
end, the membrane again was placed above the table tennis ball, supported as
before. The center of the ball was on axis of the tube. As the tube was
lowered, the membrane stretched and slipped on the ball, eventually becoming
highly stretched, as shown in Fig. 2(a). The tube was then slowly raised, with
the membrane contracting and slipping on the ball, coming to rest and forming
the invagination shown in Fig. 1.

(a)  (b)

Figure 2: (a) Forming the invagination by pushing down the PVC tube with
membrane held in place at the bottom. (b) Removing the sphere from the
invagination by pulling the rod attached to the ball.

(a) (b)

Figure 3: Experimental results. (a) Deformation $\theta$ vs. $\bar{r}$. The
inset shows the circles inscribed on the deformed membrane. (b) Principal
stretches $\lambda_{\theta}$ and $\lambda_{\phi}$ vs. $\bar{r}$. Here both the
deformation $\theta$ and undeformed distance $\bar{r}$ are normalized by the
radius of the sphere $R$.

###  2.3 Releasing the sphere

We were interested in determining the force required to release the sphere
from the invagination. We attached a solid rod (golf tee) to the ball held by
the membrane attached to PVC tube as shown in Fig. 2(b), and via pulleys and
weights, applied an increasing upward force to the rod. In this instance, it
required a force of $2.3N$ to remove the table tennis ball from the
invagination.

###  2.4 Characterizing the deformation

To characterize the deformation of the membrane, using a template, we
inscribed circles, separated by $0.25$ radians along the meridional direction,
on the stretched membrane, as shown in the inset of Fig. 3(a). The membrane
was removed from the ball, and the image on the undeformed membrane was
analyzed to characterize the deformation. The displacement $s$ of a point at
distance $r$ from the center of the undeformed membrane was measured and both
are normalized by the radius of the sphere, specifying the deformation
$\theta(\bar{r})$, as shown in Fig. 3(a). Fig. 3(b) shows the principal
stretches, as defined by Eqs. (9) and (10).

Although the two methods of realizing the invagination differ in that method 1
imposes a strong pre-stretch while method 2 does not, our measurements
indicate that the final configurations are nonetheless essentially the same.

##  3 Modeling the invagination

In this section, we present the mathematical models which describe the
deformation of the elastic membrane on the sphere involving Coulomb friction.
When the membrane loses contact with the sphere, it is still in tension and is
strictly under the sphere without making contact with it. The tension in the
membrane rapidly relaxes to zero, below which point the membrane is stress
free. In this study, we focus only on the portion of the membrane which is in
contact with the sphere. We first describe the force balance equation on the
elastic membrane, and then present the hyper-elastic models used to calculate
the forces, and finally we present our algorithm for the solution, followed by
numerical results.

###  3.1 Force balance

Our goal here is to identify the forces at play and to examine the conditions
for the stability of the invagination. The system under consideration is a
rigid sphere, partially covered by an elastic membrane as shown in Fig. 1. The
configuration there is stable; the sphere is entrapped by the membrane and
external forces must be applied to separate the two.

Figure 4: Schematic of the invagination. The force shown is exerted while the
membrane is being stretched. Only the portion of the membrane which is in
contact with the sphere is shown.

A schematic is shown in Fig. 4. The sphere with radius $R$ is at rest with its
center at the origin. Its surface exerts not only a normal force, but also a
tangential friction force, with friction coefficient $\mu$. The initially
undeformed elastic membrane is flat, regarded here for simplicity as circular
in shape, with thickness $h_{0}$. It is then stretched and placed on the
sphere so as to cover a portion of the sphere, as shown in Fig. 4. The
experimental details for producing the stable invagination are given in
Section 2.

We assume that the system has cylindrical symmetry about the $y$-axis. We
denote the points on top and bottom of the sphere as the north and south
poles, respectively; great circles passing through the poles form the
meridians along which latitude is measured, perpendicular to these are the
lines of latitude, along which the azimuth is measured. Using the standard
spherical coordinate notation, distance $s$ from the north pole is $R$
multiplied by the polar angle $\theta$ along a meridian, and points along the
azimuth are identified by the angle $\phi$. The undeformed elastic membrane is
a flat disc, where distance from the center in the plane is measured by $r$,
distance perpendicular to the plane is measured by $z$.

When the deformed membrane forming the invagination is on the sphere and the
system is in equilibrium, the net force on any element of the membrane must
vanish. There will be stress $\bm{\sigma}$ in the membrane, with principal
stress vectors in the tangential meridional, the tangential azimuthal and the
normal directions. The net force on any element of the membrane with volume
$\Delta V$ exerted by the surrounding membrane is

$\mathbf{F}_{sm}=-\int_{\Delta V}\nabla\cdot\bm{\sigma}dV,$ (2)

with components in all three principal directions. In addition to
$\mathbf{F}_{sm}$, there will be a normal force $\mathbf{F}_{N}$ per area
exerted by the sphere on the membrane, normal to the sphere, and a friction
force $\mathbf{F}_{f}$ per area in the meridional direction in equilibrium.
Due to cylindrical symmetry, there will be no friction force in the azimuthal
direction. If the stress vectors in the meridional and azimuthal directions
are multiplied by the thickness of the membrane divided by the sphere radius,
they give $T_{\theta}$ and $T_{\phi}$, the components of the tangential force
per area due to stress on the membrane in the meridional and azimuthal
directions, respectively.

The net force in the azimuthal direction is zero due to symmetry and the net
force in the normal direction is zero by Newton’s third law. Force balance in
the meridional direction requires that

$F_{m}=-{F}_{f},$ (3)

where ${F}_{f}$ denotes the (areal) friction force density, and

$F_{m}=\frac{\partial
T_{\theta}}{\partial\theta}+(T_{\theta}-T_{\phi})\cot\theta$ (4)

denotes the elastic force density along the meridian.

Friction plays a key role in the system. During the formation of the
invagination, initially, as the membrane is pulled downward, it slips on the
sphere with $s(r)$ increasing in the northern hemisphere for both methods 1
and 2. After the stretching, as the tension in the membrane is reduced, the
membrane slips again - but now towards the north pole with $s(r)$ decreasing
in some - but not all - regions on the sphere. As the membrane slips, a
tangential friction force density in the meridional direction
${F}_{f}=\mu_{k}F_{N}$ is exerted by the sphere on the membrane - the sign of
the force depends on the direction of the slip; the magnitude is determined,
in the simplest approximation, by the kinetic friction coefficient. The
membrane moves as long as the elastic force density $F_{m}$ is greater than
${F}_{f}$. As the membrane slips, $F_{m}$ decreases, and when $F_{m}<$
${F}_{f}$, the membrane comes to rest at that point. To move the membrane at
rest requires a force greater than the maximum provided by static friction.

Experimental observations indicate that the sliding membrane does not come to
rest all at once; regions near the equator continue slipping, while the region
near the north pole is at rest. Although the dynamics of the boundary between
the stationary and moving parts of the membrane is intriguing, due to the
complexity of dynamic friction [5], we will not consider it here. For
simplicity, we ignore the difference between static and kinetic friction,
assume that $\mu_{k}\approx\mu_{s}=\mu$, but stress that in the static case,
the relation ${F}_{f}=\mu F_{N}$ denotes only the maximum possible friction
force density that the sphere can provide. Explicitly, the friction force
density $F_{f}$ satisfies

$|F_{f}|\leq\mu F_{N},$ (5)

where, in our case, the normal force density is given by

$F_{N}=T_{\theta}+T_{\phi}.$ (6)

We note that $T_{\theta}$ and $T_{\phi}$ contribute to both the net normal and
net meridional force density. The force balance equation that the system must
satisfy, in the case of a stable invagination, becomes

$-\mu(T_{\theta}+T_{\phi})\leq\frac{\partial
T_{\theta}}{\partial\theta}+(T_{\theta}-T_{\phi})\cot\theta\leq\mu(T_{\theta}+T_{\phi}).$
(7)

Finally, it is interesting to consider the net force exerted by the sphere on
the entire membrane. The horizontal component of the net force is zero due to
the cylindrical symmetry, and the vertical component is

$F_{v}^{mem}=2\pi R^{2}\int_{0}^{\theta_{max}}(F_{N}\cos\theta-
F_{f}\sin\theta)\sin\theta d\theta.$ (8)

If the total vertical force on the membrane is positive, the membrane will
accelerate away from the sphere, and no stable invagination is realized. If
the total vertical force on the membrane is negative, the center or mass of
the membrane will accelerate towards the south pole. Therefore, the stable
invagination is attained only if the total vertical force $F_{v}^{mem}$is
zero.

###  3.2 Hyperelastic strain energies and forces

The components $T_{\theta}$ and $T_{\phi}$ are not independent, but are
related through stretches in the membrane. Details are provided by models of
elasticity. In order to gain a thorough understanding of the deformation
associated with the invagination, we consider three elasticity models: Mooney-
Rivlin hyperelastic model in the Hookean limit, the Neo-Hookean limit, and the
Gent model. We use the proper Cauchy stress throughout.

In spherical coordinates, the principal stretches tangential to the sphere are

$\lambda_{\theta}=\frac{\partial s}{\partial r},$ (9)

and

$\lambda_{\phi}=\frac{\sin(s/R)}{r/R},$ (10)

where $s$ is the distance to the north pole after the deformation, and $r$ is
the distance to the north pole in the underformed state. We assume that the
deformation is isochoric; then

$\lambda_{z}=\frac{1}{\lambda_{\theta}\lambda_{\phi}}.$ (11)

With the stretches $\lambda_{\theta}$ and $\lambda_{\phi}$ defined, we can
find the corresponding stresses from the various models of elasticity.

The hyperelastic incompressible Mooney-Rivlin strain energy density has the
form [6, 7, 8]

$W_{R}=C_{1}(I_{1}-3)+C_{2}(I_{2}-3),$ (12)

where

$2C_{1}+2C_{2}=G,$ (13)

and $G$ is the shear modulus, and

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\lambda_{\theta}^{2}+\lambda_{\phi}^{2}+\frac{1}{\lambda_{\theta}^{2}\lambda_{\phi}^{2}},$
(14) $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\lambda_{\theta}^{2}\lambda_{\phi}^{2}+\frac{1}{\lambda_{\theta}^{2}}+\frac{1}{\lambda_{\phi}^{2}}$
(15)

are the two invariants of the left Cauchy-Green deformation tensor. If
$C_{2}=0$, it gives the Neo-Hookean energy density

$W_{NH}=\frac{1}{2}G(\lambda_{\theta}^{2}+\lambda_{\phi}^{2}+\frac{1}{\lambda_{\theta}^{2}\lambda_{\phi}^{2}}-3).$
(16)

In the limit of small strains, we have the Hookean energy density,

$W_{H}=2G((\lambda_{\theta}-1)^{2}+(\lambda_{\phi}-1)^{2}+(\lambda_{\theta}-1)(\lambda_{\phi}-1)).$
(17)

The Gent elasticity model takes into account the finite extensibility [9],
with

$W_{G}=-\frac{GJ_{m}}{2}\ln(1-\frac{I_{1}-3}{J_{m}}),$ (18)

where $I_{m}=J_{m}+3$ represents the maximum value of $I_{1}$ allowed by the
system in the absence of failure. In the limit as $I_{m}\rightarrow\infty$,
the Gent model reduces to the Neo-Hookean model.

The tangential force density components $T_{\theta}$ and $T_{\phi}$, in terms
of the stretches, are

$T_{\theta}=\frac{h_{0}}{R}\frac{1}{\lambda_{\phi}}\frac{\partial
W}{\partial\lambda_{\theta}},$ (19)

and

$T_{\phi}=\frac{h_{0}}{R}\frac{1}{\lambda_{\theta}}\frac{\partial
W}{\partial\lambda_{\theta}}.$ (20)

Substitution of Eqs. (19) and (20) into Eq. (7) gives the inequality in terms
of $s(r)$ whose solution provides the detailed description of the elastic
deformation. We note that the shear modulus plays no role in the invagination.

###  3.3 Numerical algorithm and results

Our goal is to find the stable deformation of the membrane $s(r)$, for which
Eq. (7) must be satisfied. This requires a non-standard approach, which we now
describe. All the lengths are nondimensionalized by the radius of the sphere
$R$, so $s/R=\theta$, and $r/R=\bar{r}$. We first prescribe the portion of the
membrane in contact with the sphere by $\bar{r}_{\max}$, and the point on the
sphere $\theta_{\max}$, where the membrane loses contact. In this example, we
have set $\bar{r}_{\max}=0.9$, $\theta_{\max}=1.8$, which are in rough
agreement with experiments. We also prescribe the initial displacement
$\theta_{0}(\bar{r})$ of the membrane. In particular, we use an initial
displacement of the membrane which is in the neighborhood of the final
displacement, which has been suggested by experiments. Given a friction
coefficient $\mu$, we implement the following psudeodynamics

$\theta_{i}^{n+1}=\theta_{i}^{n}+\Delta\theta(\bar{r}_{i})dt,$ (21)

where $\theta_{i}^{n}$ is the position of membrane at undeformed position
$\bar{r}_{i}$ at time step $n$, and

$\Delta\theta(\bar{r}_{i})=\left\\{\begin{tabular}[]{ll}$0,$&if
$|F_{m}|\leq\mu F_{N},$\\\ $F_{m}-\mu F_{N},$&if $F_{m}>\mu F_{N},$\\\
$F_{m}+\mu F_{N},$&if $F_{m}<-\mu F_{N}.$\end{tabular}\right.$ (22)

That is, the point on the membrane does not move if $F_{m}\,$is less than the
maximum frictional force, and the point moves either away from or towards the
north pole if $F_{m}$ overcomes the maximum frictional force, and the
direction of motion depends on the sign of $F_{m}$. We have used 100 equally
spaced points $\bar{r}_{i}$ between $0$ and $\bar{r}_{\max}$. All derivatives
are approximated with a centered difference scheme.

After equilibrium is reached, the total vertical force exerted on the membrane
by the sphere can be evaluated by Eq. (8). If $F_{v}^{mem}$ is different from
zero, we adjust the value of the friction coefficient $\mu$. Explicitly, we
increase $\mu$ if $F_{v}^{mem}>0$, and decrease $\mu$ otherwise, until the
total vertical force on the membrane vanishes.

We remark that the equilibrium solution is highly sensitive to the initial
condition, so is the friction coefficient which stabilizes the invagination.

Fig. 5 shows the numerical results from the Hookean, Neo-Hookean and Gent
models. The results from three different models are qualitatively the same,
with the main difference that they require different values of friction
coefficient to stabilize the invagination. Among the three models, the Hookean
requires the largest friction coefficient, and Neo-Hookean requires the
smallest friction coefficient. The boundary of the two background colors on
the figures marks the change of direction of the tangential force from
elasticity, or equivalently the change of direction of the frictional force.
The boundary also coincides with the location where the curve
$\theta(\bar{r})$ describing the deformation curve changes its concavity (left
column of Fig. 5), where the two stretches reach their maxima(middle column of
Fig. 5). We note that it is not at the north pole where the membrane is
maximally stretched, but rather at an angle $\theta$ close to $\pi/4$, where
the membrane pulls the neighboring regions both sides towards each other. The
sphere pushes the membrane upward with increasing force from the north pole to
the point with maximum stretch, where the force changes direction and pulls
the membrane downward towards the point where the membrane loses contact with
the sphere, as shown in right column of Fig. 5. Remarkably, in the equilibrium
configuration, the tangential force is in balance with the maximum friction
force everywhere, except at a finite number of discrete points where the
stretches reach a local extremum and the friction force changes direction. Any
decrease of the friction coefficient $\mu$ from its thus found equilibrium
value will render the invagination unstable.

Hookean

NeoHookean

Gent

Figure 5: Top row: results from Hookean model with $\mu=1.47$. Middle row:
results from neo-Hookean model with $\mu=0.99$. Bottom row: results from Gent
model with $J_{m}=50$, $\mu=1.1$. Left column: deformation $\theta$ vs.
$\bar{r}$. Middle column: principal stretches $\lambda_{\theta}$ and
$\lambda_{\phi}$ vs. $\bar{r}$. Right column: elastic force density along
meridional direction $F_{m}$ vs. $\bar{r}$.

##  4 The situation in 2D

It is interesting to ask whether such an invagination can be realized with an
elastic membrane and a circular cylinder. A quick experiment with a membrane
on a circular cylinder shows that the answer is no. If the cylinder’s axis is
horizontal, the membrane can’t stay on the cylinder under the midplane without
an external force, which is in stark contrast with the spherical case, where
the hoop stress from azimuthal direction can keep the membrane against the
sphere below the equator. Second, in cylindrical geometry, the force balance
equation reads

$-\mu T_{\theta}\leq\frac{dT_{\theta}}{d\theta}=-F_{f}\leq\mu T_{\theta},$
(23)

where $T_{\theta}$ is the force per length. The solution for $T_{\theta}$ at
any angle $\theta$ is bounded by

$T_{0}\exp(-\mu\theta)\leq T_{\theta}\leq T_{0}\exp(\mu\theta),$ (24)

where $T_{0}$ is the tension at the top of the cylinder where $\theta=0$. The
total vertical force on the membrane from the cylinder can be written as

$\displaystyle F_{v}^{mem}$ $\displaystyle=$ $\displaystyle
R\int_{0}^{\theta_{max}}(T_{\theta}\cos\theta-F_{f}\sin\theta)d\theta$ (25)
$\displaystyle=$ $\displaystyle
R\int_{0}^{\theta_{max}}(T_{\theta}\cos\theta+\frac{dT_{\theta}}{d\theta}\sin\theta)d\theta$
$\displaystyle=$ $\displaystyle
R\int_{0}^{\theta_{max}}d(T_{\theta}\sin\theta)=RT_{\theta}\sin\theta_{\max}.$

From Eq. (24), since $\theta_{\max}\in(0,\pi)$, the force on the membrane will
be always positive. Thus there is no elastically stabilized invagination in
2D, where the hoop stress from azimuthal direction is missing.

##  5 Discussion and Summary

By conducting experiments using regulation size table tennis balls and dental
dams, we have realized elastically stabilized spherical invaginations. We used
two approaches to form the invagination; they gave similar deformations of the
elastic membrane. Friction plays a key role in stabilizing the system.

Using $|F_{f}|\leq\mu F_{N}$, we obtained the governing inequality for the
elastic deformation of the membrane on the sphere. The forces at play are
expressed in terms of stretches in the membrane. We used three elastic models
to describe the invagination; crude agreement with experiments was found.
Although the inequality was used to evolve the dynamics, remarkably, the
solution obeys the equality at almost every point.

Both the experimental data and the numerical solutions indicate that the
maximum stretch of the membrane occurs at $\theta\simeq\pi/4$, where the both
frictional force and vertical component of the vertical force on the membrane
change direction. Due to the highly nonlinear nature of the problem, the
equilibrium solution and the required minimum friction coefficient are
strongly dependent on initial conditions, or equivalently, on the process
establishing the invagination.

We note the close similarity of the invagination of the membrane and the table
tennis ball and the invagination of the membrane and a coin. In the latter
case, the membrane was nearly uniformly stretched on top of the coin, and the
stretch slowly decreases to a relaxed state under the coin. The invagination
was stabilized by the friction force between the vertical edge of the coin and
the membrane. Here, the vertical downward force on the membrane is from the
lower part of the sphere, below the latitude where the membrane is maximally
stretched.

When releasing the ball from the membrane, there is no vertical force to
overcome, since the the total vertical force on the sphere exerted by the
membrane is zero. However, there is hoop stress must be overcome. The
retracting force must be sufficiently large to deform the membrane and
increase the radius of the bottom opening to the radius of the sphere. The
sphere may also be released if the elastic and frictional properties of the
membrane changed. In addition to the experiments described here, we have also
realized spherical invaginations using other membranes (balloons, condoms,
liquid crystal elastomers) and spherical rigid bodies (marbles, push pins,
ball chains).

Here we have only focused on the basic mechanism of the invagination and
considered only the portion of the membrane in contact with the sphere.
Modeling the entire membrane would give a more complete description of the
invagination.

## Acknowledgments

This work was supported by the Office of Naval Research through the MURI on
Photomechanical Material Systems (ONR N00014-18-1-2624) and Air Force contract
FA8649-20-C-0011 as part of the STTR AF18B-T003 Electronically Dimmable Eye
Protection Devices (EDEPD) program.

## References

  * [1] Stillwell, W. An Introduction to Biological Membranes. Elsevier Scientific (2016).
  * [2] Meng, F., Doi, M., Ouyang, Z. et al. The ‘Coin-Through-the-rubber’ Trick: An Elastically Stabilized Invagination. J Elast. 123, 43–57 (2016).
  * [3] L. Euler, Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes, sive solutio problematis isoperimetrici latissimo sensu accepti, Additamentum I (De curvis elasticis) 245–320 (1744).
  * [4] Varenberg M., Varenberg A. Table Tennis membrane: Tribological Characterization, Tribol. Lett. 47, 52-56 (2012).
  * [5] Liu, Y.F. et al. Experimental comparison of five friction models on the same test-bed of the micro stick-slip motion system, Mech. Sci. 6, 15-28 (2015).
  * [6] Mooney, M., A theory of large elastic deformation, Journal of Applied Physics, 11(9), 582–592 (1940).
  * [7] Rivlin, R. S., Large elastic deformations of isotropic materials. IV. Further developments of the general theory, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 241(835), 379–397 (1948)
  * [8] Audoly, B., Pomeau, Y.: Elasticity and Geometry: From Hair Curls to the Non-linear Response of Shells. Oxford University Press, Oxford (2010).
  * [9] Gent, A.N., A new constitutive relation for membrane, Rubber Chemistry and Technology, 69, 59–61 (1996).