**2. Materials and Methods**

We utilized de-identified and aggregated large-scale data developed by Facebook Data for Good platform to identify population mobility trends during the COVID-19 crisis [33]. These data are called Movement Data and they aggregate information from people using the Facebook app on their mobile phones with location history turned on to show movements between two points. The computation system of this data acquires the most common location of the Facebook app user using Bing tile map Level 13 (e.g., approximately 4.9 × 4.9 km) [34] within the first time window and the most common location in the second time window for each day. The centroids of the starting and ending Bing tiles are assigned to the person's movement vector. Then, these vectors are aggregated into the administration level 4 boundary (e.g., township), which is spatially equivalent to MCDs in Maryland and CCDs in California. MCDs and CCDs are administrative county subdivisions for which the Census Bureau establish and provide subcounty statistics [35]. CCDs are equivalent geographic entities to MCDs in US states where MCDs do not exist or have been unsatisfactory for comparing statistical data [35]. The data characterized mobility trend as the percent change in the observed number of users traveling into the administrative area (MCD or CCD) for the same time window and the same day of the week compared to the baseline period at MCD (CCD) level so the mobility data indicate movement of users across regions at this spatial level. The baseline number of users moving into an MCD (CCD) was calculated as the average number of users moving for the same daytime window and day of the week between 26 February 2020 to 10 March 2020. The percent change in the number of users between a given day and its baseline is calculated as follows [33]: percent change = (*c<sup>t</sup>* − µ*baseline*,*t*/ µ*baseline*,*<sup>t</sup>* + ǫ , where *c* is the number of users moving into a MCD (CCD) for day *t*, µ*baseline*,*<sup>t</sup>* is the mean of the number of users traveling into the same MCD over the same time interval on the same day of the week of day *t* during the baseline period, and ǫ is a small value, in this case 1. The distance traveled by users for each MCD (CCD) was calculated as the distance of the movement vectors linking the centroids of the starting and end Bing tiles of all users who traveled into that MCD (CCD). We used the daytime window (8 a.m.–4 p.m.) for our analysis.

The Movement Data of our study regions were linked to the geographic information systems (GIS) data of MCD in Maryland and CCD in California provided by the US Census Bureau. We analyzed the 230 administrative areas (MCDs) out of 290 total administrative areas in the state of Maryland, US, for which the Movement Data were available, in order to examine mobility trends for the period from 10 March to 24 April 2020. For California, we analyzed 341 CCDs out of 397 CCDs, for which the Movement Data were available. Percent change in the number of users moving was not observed for some MCDs when MCDs are smaller than a Bing tile so the user's location at MCD and CCD level cannot be specified. Data are not provided for county subdivisions where the number of observed users is smaller than 10 users to protect users' privacy. County subdivisions with no observation for users' moving between pairs of county subdivisions for 70% or more of the days in the study period were excluded in the analysis. As a result, 76 MCDs in Maryland and 241 CCDs were included in our main analysis.

Vegetation level to indicate green space was estimated by the Enhanced Vegetation Index (EVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) product MOD13Q1, which is a 16-day composite image at 250-meter resolution. The EVI is an advanced version of the Normalized Different Vegetation Index (NDVI). The NDVI is calculated as near-infrared radiation minus visible radiation divided by near-infrared radiation plus visible radiation (i.e., NDVI = (NIR − RED)/(NIR + RED)). The index ranges from −1 to +1 with higher values indicating denser vegetation and −1 indicating waterbody features (NASA, 2018). While the EVI's calculation is similar to NDVI, it corrects for some distortions in the reflected light caused by the particles in the air, ground cover below the

vegetation, and the saturating effects of areas with large amount of chlorophyll such as rainforests [36]. We calculated the average EVI for each study region (MCD, CCD) using the EVI pixel values within and surrounding the MCD and CCD boundary, for 1 January 2020 through 21 March 2020 to represent vegetation level in the study regions.

The MODIS Land Cover Type Product (MCD12Q1) [37] was used to estimate urbanicity of each study region. The number of 'Urban and Built-up Lands' pixels based on the University of Maryland legend and class definition was divided by the total number of pixels within the MCD or CCD boundary was calculated as the percent impervious area.

We obtained the park data in California from Esri Data & Maps [38]. These data include parks and forests at national, state, county, and local levels (e.g., city-scale parks, pocket parks, playgrounds, etc.) [38]. We obtained the Maryland State Parks data provided by the OpenStreetMap [39]. The OpenStreeMap is a global dataset including open user-generated street maps, geographical features, and built environment, which have been used for a wide range of studies [40]. This Maryland State Parks data, developed by the Baltimore County Government, is a GIS shape file and includes geographical polygons for several types of green space: state and national parks or forest, hiking trails, natural resource management areas for recreational activities (e.g., fishing, hunting, wild animal observation, etc.) and preservation of environmental resources, and wildlife management areas [41]. Data of parks at local levels in Maryland were obtained from Esri Data & Maps [38]. Hereafter, we use the term parks to refer to all these types of areas. The continuous variable for parks did not have a normal distribution. Thus, we considered a categorical variable of presence of parks within MCDs and CCDs (i.e., MCDs with state parks vs. MCDs without state parks within their boundary).

We obtained the GIS file of food retail establishments and hospitals from various sources. We obtained the data for food stores (2017–2018) and restaurants (2019) in Maryland from the Maryland Food System Map (JHSPH) developed by the Johns Hopkins Center for a Livable Future [42]. The food stores data included attributes for grocery stores, supermarkets, gas stations, and pharmacies. The GIS data of hospitals (acute, general, and special) licensed by the Maryland Department of Health and Mental Hygiene Office of Health Care Quality were obtained from Maryland's Mapping & GIS Portal [43]. Using these datasets, we calculated the sum of the number of food stores (grocery stores, supermarket, gas station, pharmacy), restaurants, pharmacies, and hospitals for each study MCD in Maryland. For California, the data on hospitals (2020) were obtained from Esri Data & Maps [38] and the data for pharmacies (2019) were obtained from OpenStreetMap [39]. The data for food retail establishments were not available for California in this study due to the lack of data for many CCDs across California.

We calculated statistics such as the first quartile (Q1), third quartile (Q3), mean, and minimum values of the daily percent changes in number of people moving between pairs of study MCDs (or CCDs) between 31 March and 24 April 2020 (i.e., after stay-at-home order) for Maryland and between 31 March and 19 April for California to characterize the mobility trends. Using linear regression analysis, we analyzed the relationships between the local vegetation level (i.e., EVI), presence of parks, and percent changes in mobility. A linear regression analysis was used for these statistics and parks, EVI, and urbanicity to examine if mobility patterns during the early stage of COVID-19 pandemic differed by green space (i.e., incorporating state parks or EVI level). We applied several different statistical models with different sets of confounders. The Q3 of mobility changes was used as a dependent variable as it represented the best normal distribution in Maryland (Supplementary Figure S1), while it was slightly skewed in California. We examined if results for green space and mobility trends were confounded by urbanicity (e.g., population, percent impervious area). We conducted regression analyses separately for each state.
