Estimation of Spatial Distribution and Fluid Fraction of a Potential Supercritical Geothermal Reservoir by Magnetotelluric Data: A Case Study From Yuzawa Geothermal Field, NE Japan

Fluids within the Earth's crust may exist under supercritical conditions (i.e., >374°C and >22.1 MPa for pure water). Supercritical geothermal reservoirs at depths of 2–10 km below the surface in northeastern (NE) Japan mainly consist of magmatic fluids that exsolved from the melt during the course of fractional crystallization. Supercritical geothermal reservoirs have received attention as next‐generation geothermal resources because they can offer significantly more energy than that obtained from conventional geothermal reservoirs found at temperatures <350°C. However, the spatial distribution and fluid fraction of supercritical geothermal reservoirs, which are required for their resource assessment, are poorly understood. Here, the magnetotelluric (MT) method with electrical resistivity imaging is used in the Yuzawa geothermal field, NE Japan, to collect data on the fluid fraction and spatial distribution of a supercritical geothermal reservoir. The collected MT data reveal a potential supercritical geothermal reservoir (>400°C) with dimensions of 3 km (width) × 5 km (length) at a depth of 2.5–6.0 km below the surface. The estimated fluid fraction of the reservoir is 0.1%–4.2% with salinity values of 5–10 wt%. The melt is also imaged below the reservoir, and based on the resistivity model; we develop a mechanism for the evolution of the supercritical geothermal reservoir, wherein upwelling supercritical fluids supplied from the melt are trapped under less permeable silica sealing and accumulate there.

. Supercritical geothermal reservoirs at a depth of 2-10 km in northeastern (NE) Japan, which is the subject of this study, are considered to consist mainly of magmatic fluids (Tsuchiya et al., 2016). In fact, several deep drilling expeditions have found magmatic fluid in a supercritical state. An exploration well  in the Kakkonda geothermal area in NE Japan found magmatic fluid at 500°C at a depth of 3,729 m Ikeuchi et al., 1998), and the fluid was in a supercritical state (Saishu et al., 2014;Tsuchiya & Hirano, 2007). Similarly, an exploratory well drilled by the Iceland Deep Drilling Project (IDDP) also encountered a supercritical fluid (450°C) in the Krafla volcanic system, Iceland (Elders et al., 2014;Scott et al., 2015).
Supercritical geothermal reservoirs have recently gained attention as next-generation geothermal resources (Friðleifsson et al., 2020;Okamoto et al., 2019;Parisio et al., 2019;Reinsch et al., 2017;Stimac et al., 2017;N. Watanabe et al., 2017). Various countries, including Japan, the USA, New Zealand, and Italy, are beginning to consider the development of power plants using supercritical fluids (Reinsch et al., 2017). This newfound attention comes from the fact that supercritical fluids have significant advantages as a source of energy: a higher enthalpy per unit mass and a lower fluid viscosity than conventional geothermal fluids with temperatures <350°C (Reinsch et al., 2017). Power generation systems using supercritical fluids are therefore able to achieve a large output from a single power plant. A test of the capacity of supercritical fluids (wellhead temperature of 450°C) from IDDP-1 in Krafla (Iceland) revealed that the energy output was 10-fold higher than that of conventional wells (Elders et al., 2014).
Gaining an in-depth understanding of the spatial distribution and fluid fraction of supercritical geothermal reservoirs is necessary for their assessment as a next-generation energy resource (Reinsch et al., 2017;Stimac et al., 2017). Deep drilling can be used to estimate the number of available resources of supercritical fluids Elders et al., 2014;Fournier, 1999;Friðleifsson et al., 2020;Ikeuchi et al., 1998). However, as a method to assess supercritical geothermal reservoirs, deep drillings are expensive (Elders et al., 2014). Owing to the limited number of deep drilling points, the distribution and fluid fraction of supercritical geothermal reservoirs remain poorly understood. Therefore, geophysical methods (e.g., electromagnetic and seismic methods) are required to gain insights into the spatial distribution and fluid fraction of supercritical geothermal reservoirs (Piana Agostinetti et al., 2017;Reinsch et al., 2017).
In this study, we use the magnetotelluric (MT) method with electrical resistivity imaging to estimate the spatial distribution and fluid fraction of a supercritical geothermal reservoir in NE Japan. The MT method is suitable for measuring the spatial distribution and fluid fraction of subsurface fluids because electrical resistivity is sensitive to the existence of fluids (Comeau et al., 2020;Ogawa et al., 2014;Unsworth et al., 2005;Wannamaker et al., 2009). The fluid fraction of a supercritical geothermal reservoir can be obtained from a resistivity model using MT data. For the estimation of the fluid fraction, we consider a two-phase mixing law (Hashin & Shtrikman, 1962), assuming that the supercritical geothermal reservoir consists of solid rock and supercritical fluid in the solid-rock pore space and that the supercritical fluid saturates the pore space. The resistivity obtained using the MT method is bulk resistivity, which includes contributions from the solid rock and supercritical fluid in the pore space (Chave & Jones, 2012;Hashin & Shtrikman, 1962). We can obtain the fluid fraction of a supercritical geothermal reservoir from the bulk resistivity using the mixing law if the resistivities of a supercritical fluid and a solid rock are known.
We consider magmatic supercritical fluids as NaCl-H 2 O fluids because magmatic fluids mainly consist of water and chloride salts (Blundy et al., 2021;Heinrich, 2005;Monecke et al., 2018;Richards, 2011;Sillitoe, 2010). The electrical resistivity of NaCl-H 2 O fluids has been studied extensively in the relationship between subsurface resistivity models and geological structures (Bannard, 1975;Nono et al., 2020;Quist & Marshall, 1968;Sakuma & Ichiki, 2016;Sinmyo & Keppler, 2016). Bannard (1975) reported that the electrical conductivity (the reciprocal quantity of resistivity) of NaCl-H 2 O fluids strongly depends on NaCl concentrations, temperature, and pressure based on measurements of NaCl-H 2 O fluids for a wide range of pressure, temperature, and NaCl conditions, including supercritical conditions (up to 200 MPa, 525°C, and 25 wt% NaCl concentration). The measurements showed that the conductivity of the NaCl-H 2 O fluid increases with increasing NaCl concentration under a fixed temperature and pressure condition (e.g., 0.6 S/m and 115 S/m for 0.06 wt% and 24.6 wt% NaCl concentrations, respectively at 25 MPa and 290°C). The measurements also revealed that the conductivity of NaCl-H 2 O fluids under a fixed NaCl concentration and pressure condition increases with increasing temperature between 0°C and 300°C, reaches a maximum near 300°C, and decreases as temperature increases between 300°C and 525°C 10.1029/2021JB022911 3 of 21 (Bannard, 1975). This result indicates that supercritical fluids tend to be more electrically resistive compared to fluids at 300°C. Recent studies have also confirmed the decrease in conductivity in the supercritical region (Kummerow et al., 2020;Nono et al., 2020). N. Watanabe et al. (2021) developed a viscosity-dependent empirical model to calculate the conductivity of NaCl-H 2 O fluids based on Bannard (1975) measurements. We use this viscosity-dependent empirical model to calculate the conductivity of supercritical fluids. The solid rock in supercritical geothermal reservoirs in NE Japan is considered to be granite Pastor-Galán et al., 2021). We use the resistivity of dry granite reported by Kariya and Shankland (1983) for the resistivity of the solid rock.
NE Japan can be classified as a typical subduction zone, where the Pacific plate subducts beneath the land area at a rate of ∼10 cm/year (Hasegawa et al., 1991). The dehydration of the subducting slab in the mantle wedge results in the upward migration of fluids (Tatsumi, 1989). Since the presence of water lowers the melting temperature of rocks, the upwelling hot fluids can result in partial melting of the mantle wedge and the crust (Iwamori, 1998). The partial melting in NE Japan was detected by low-velocity anomalies of S wave (so-called "hot fingers") and these hot fingers are distributed over NE Japan at intervals of several tens of kilometers (Tamura et al., 2002). Magmatic fluids exsolved from the melt are considered to contribute to the formation of supercritical geothermal reservoirs at a depth of 2-10 km in NE Japan (Tsuchiya et al., 2016). Therefore, the existence of these hot fingers indicates that a large number of supercritical geothermal reservoirs may exist in NE Japan (Okamoto et al., 2019).
We select the Yuzawa geothermal field in NE Japan as the target area. Hot fingers exist around this geothermal field (Tamura et al., 2002). Other geophysical and geochemical data (e.g., thermoluminescence, borehole temperature, and seismic data) also suggest that supercritical fluids may exist at a depth of a few kilometers (New Energy and Industrial Technology Development Organization, 1990;Nunohara et al., 2021;Okada et al., 2014). We conducted an MT survey to reveal the spatial distribution and fluid fraction of a potential supercritical geothermal reservoir in this field. Any information obtained on a supercritical geothermal reservoir in this field is likely to be useful for the estimation of spatial distribution and fluid fraction of other supercritical geothermal reservoirs in NE Japan (e.g., Naruko, Onikobe, and Kakkonda areas). Here, we first provide a brief introduction of the Yuzawa geothermal field, MT method, and a two-phase mixing law. Then, we present our findings on the estimation of the spatial distribution of a potential supercritical reservoir and its fluid fraction.

Yuzawa Geothermal Field in NE Japan
The Yuzawa geothermal field in NE Japan contains Wasabizawa geothermal power plant, the fourth-largest geothermal power plant in Japan (Nunohara et al., 2021). This field is located in the inner part of the Sanzugawa caldera (Figures 1a and 1b). The Sanzugawa caldera, which was first formed ca. 3 Myr ago, has a size of 30 km north-south and 20 km east-west (Takeno, 2000). A Quaternary volcano, the Takamatsu volcano, or Takamatsu-dake, is found in this geothermal field (red triangle in Figures 1b and 1c). The Takamatsu volcano is located 40 km west of the volcanic front (Tamanyu et al., 1998). The volcano is composed of calc-alkaline andesite to dacite. The age of this volcano was estimated to be 0.2-0.3 Myr by thermo-luminescence and K-Ar methods (Umeda et al., 1999). The magma below the volcano is considered to be the heat source of this geothermal field (Takeno, 2000). Geothermal features, such as hot springs and hydrothermally altered rocks, have been observed extensively at the surface in this geothermal field, indicative of a well-developed hydrothermal system (Nunohara et al., 2021). Two geothermal power plants (Uenotai: 28,500 kW and Wasabizawa: 46,200 kW) are currently in operation using a hydrothermal system at temperatures <300°C (magenta stars in Figure 1c).
Seismic imaging of NE Japan suggests that upwelling fluids exist beneath this geothermal field at a depth of up to several kilometers (Nakajima et al., 2001;Okada et al., 2014;Tamura et al., 2002). The high He 3 /He 4 in this field also indicates that fluids upwell to the surface from the upper mantle (Horiguchi et al., 2010;Kita et al., 1992). The temperature profile of a well, N63-MS-6 (blue square in Figure 1c) drilled by the New Energy and Industrial Technology Development Organization (NEDO), is a conduction-dominated type with a temperature gradient of 170°C/km (New Energy and Industrial Technology Development Organization, 1990). This thermal gradient suggests that temperatures can reach >400°C at a depth of 2.5 km below ground surface (bgs). These geophysical observations and borehole temperature measurements strongly indicate the possible existence of a supercritical geothermal reservoir (>400°C) in this field.

MT Method
The MT method is a geophysical tool used to map the electrical resistivity structure of subsurfaces (Araya Vargas et al., 2019;Árnason et al., 2010;Becken et al., 2011;Bedrosian et al., 2018;Di Paolo et al., 2020;Heinson et al., 2018;Heise et al., 2008;Hyndman & Shearer, 1989;Ichihara et al., 2016;Ingham et al., 2009;Le Pape et al., 2015;Moorkamp et al., 2019;Wise & Thiel, 2020). This method uses surface measurements of naturally occurring low-frequency electromagnetic variations to investigate underground resistivity structures (Chave & Jones, 2012). Each MT station records two components of the electric field (E x and E y ) and three components of the magnetic field (H x , H y , and H z ). The x-, y-, and z-directions respectively correspond to the geographic north, east, and vertical depths. The z-direction is pointing down. The orthogonal components of the horizontal electric and magnetic fields are linearly related by the complex-valued impedance tensor Z in the frequency domain as follows: (1) The tipper T (or vertical magnetic field transfer function) relates the vertical magnetic fields (H z ) to the horizontal magnetic field components in the frequency domain with the following relation: and T is dimensionless. Z and T carry information about underground resistivity structures and therefore depend on the measurement location r (x, y, z) and frequency f in Hz (Chave & Jones, 2012). The Z can be expressed in terms of the apparent resistivity and phase for data interpretation (Chave & Jones, 2012). The Z can also be expressed as a phase tensor and phase tensor is independent of the distorting effects produced by localized near-surface resistivity heterogeneities (Caldwell et al., 2004): where X and Y are the real and imaginary parts of the Z, respectively. The phase tensor Φ can be graphically represented as an ellipse with principal (ellipse) axes, Φ max and Φ min (Caldwell et al., 2004): The Φ 2 indicates the magnitude of the phase tensor response. In a one-dimensional subsurface resistivity structure, the decrease in resistivity with increasing depth is indicated by a value of Φ 2 > 45° and the increase in resistivity with increasing depth is indicated by a value of Φ 2 < 45°.
The observed MT data can be converted into an underground resistivity model using inversion analysis (Constable et al., 1987;Farquharson, 2008;Kelbert et al., 2014;Newman, 2014;Ogawa & Uchida, 1996;Siripunvaraporn & Egbert, 2009). We use WSINV3DMT for the three-dimensional (3-D) inversion of the observed data (Siripunvaraporn & Egbert, 2009). WSINV3DMT searches for the smoothest model subject to the desired level of data misfit χ * . Mathematically, this equals to searching for the stationary point of an unconstrained functional U(m, λ): where m is the resistivity model, m 0 is the prior model, d is the observed data, F[m] is the forward response, C m is the model covariance, C d is the data covariance matrix, and λ −1 is the Lagrange multiplier. The first term in the function U is the model roughness and the second term is the data misfit weighted by the Lagrange multiplier λ −1 .
Since U is a function of both m and λ, the search for the stationary points is not straightforward (Siripunvaraporn & Egbert, 2009). Alternatively, WSINV3DMT minimizes a penalty functional W λ : This is because when λ is fixed, and obtain the same result (Siripunvaraporn & Egbert, 2009). By minimizing W λ with a series of λ, the stationary point of U can be found (i.e., λ can be found such that the data misfit is 2 * ). For further details on WSINV3DMT, see Siripunvaraporn and Egbert (2009).

Two-Phase Mixing Law
We consider a two-phase mixing law proposed by Hashin and Shtrikman (1962) to estimate the fluid fraction of a potential supercritical geothermal reservoir from an electrical resistivity model. To apply the mixing law, we assume that the supercritical geothermal reservoir consists of solid rock and supercritical fluid in the solid-rock pore space and that the supercritical fluid saturates the pore space. Hashin and Shtrikman (1962) proposed the upper and lower bounds for effective conductivity. The Hashin-Shtrikman upper bound model (HS+; Hashin & Shtrikman, 1962), which assumes complete connectivity: where 1 is the conductivity of the solid rock, 2 is the conductivity of the supercritical fluid, and 2 is the volumetric fluid fraction of the supercritical fluid (i.e., porosity). Equation 7 indicates that we can obtain the fluid fraction of a supercritical geothermal reservoir from the bulk conductivity using the mixing law if conductivity values of a supercritical fluid and a solid rock are known. With the HS+ model, the contribution of solid rock conductivity to bulk resistivity is limited if the conductivity of the solid rock is much smaller than the conductivity of the fluid (Ogawa et al., 2014). The Hashin-Shtrikman lower bound model (HS-; Hashin & Shtrikman, 1962) is presented as follows: .
The HS-model assumes that the fluid is included in isolated pockets. This study only considers the Hashin-Shtrikman upper bound model (HS+) because the HS+ has provided a reasonable estimation of a fluid fraction of crustal rocks in NE Japan (Ogawa et al., 2014). We use a viscosity-dependent empirical model developed by N. Watanabe et al. (2021) to calculate the conductivity of supercritical fluids ( 2 ). We consider that the solid rock in supercritical geothermal reservoirs in NE Japan consists of granite  and use the conductivity of dry granite reported by Kariya and Shankland (1983) for the conductivity of the solid rock ( 1 ).

Results
We collected MT data at 30 measurement sites to image resistivity structures in the Yuzawa geothermal field ( Figure 1c). Since the objective of this study is to image supercritical geothermal reservoirs at depths of several kilometers, the range of the MT measurement area is 15 km in the horizontal direction. The Φ 2 is higher than 45° at 20 Hz and lower than 45° at 2.5 Hz around Takamatsu volcano and power plants (Figures 2a and 2b), suggesting the existence of near-surface conductors and resistors below the near-surface conductors. Apparent resistivity data at MT sites around Takamatsu volcano and power plants exhibit low values <50 Ωm at a high frequency of 100 Hz, which also suggests the existence of the near-surface conductors around Takamatsu volcano and power plants (Figures S1, S2, and S4 in Supporting Information S1). The Φ 2 is higher than 45° at 0.3 Hz around the northwest of the Takamatsu volcano, which implies the presence of a deep conductor ( Figure 2c). The impedance and tipper data for each observation site are shown in Figures S1-S6 in Supporting Information S1.
Inverse modeling with WSINV3DMT can convert the observed MT data into 3-D resistivity structures (Siripunvaraporn & Egbert, 2009;Siripunvaraporn et al., 2005). The initial resistivity model for inverse modeling consists of the subsurface (100 Ωm), ocean (0.3 Ωm), and air (10 10 Ωm). We tested other initial models with homogeneous subsurface structures of 1, 10, and 1,000 Ωm. The initial model with homogeneous subsurface structures of 100 Ωm obtained the lowest root-mean-square (RMS) final misfit (the final RMS misfit = 6.37, 3.33, 2.34, and 3.55 for the initial models of 1, 10, 100, and 1,000 Ωm, respectively). The prior resistivity model was set to be the same as the initial resistivity model. The model mesh is discretized into 74 × 74 × 70 cells in the x-, y-, and z-directions, respectively, containing boundary cells. The horizontal cell size is kept constant at 300 m inside the area of interest (−10 km < x and y < 10 km) and then is increased toward the boundary. For depths between −1.3 km below sea level (bsl) and 0 km bsl (depths <0 km bsl indicate above the sea level in our definition), the vertical cell size is set to values ranging from 20 to 50 m. For depths >0 km bsl, the vertical cell size is increased as the depth increases. We performed the inverse modeling on a computer (@Xeon 3.10 GHz Gold 6254 central processing 7 of 21 unit; Intel Corp.) with 3 TB of RAM. The impedance tensor and tipper between frequencies of 0.001-100 Hz (17 frequencies) are used for inverse modeling.
Error floors for the impedance tensor are set to 5% of | | 1 ∕2 ; error floors for the tipper are set to 15% of each component. The smoothing parameters for C m are set to = 5 , = 0.1, = 0.1, and = 0.3. We selected these smoothing parameters because the inverse modeling with these smoothing parameters obtained the lowest final RMS misfit compared to other smoothing parameters. For details on the smoothing parameters of C m , see Siripunvaraporn and Egbert (2009). The RMS misfit for the initial model was 10.2. The inversion analysis obtained a resistivity model with an RMS misfit of 2.34 after three iterations ( Figure 3). The data fits are good for almost all sites and frequencies (Figures S1-S6 in Supporting Information S1).
3-D resistivity structures obtained from the inverse modeling identify near-surface conductors at a depth between −1.3 and 0.5 km bsl, as well as a deep conductor C1 at a depth between 1 and 12 km bsl (Figure 3). The resistivity of the near-surface conductors is 5-20 Ωm, which is lower than that of the surrounding rock of 100 Ωm. The near-surface conductors are widely distributed in this area. The resistivity of C1 is 5-15 Ωm. The size of C1 is 3 km (width) × 5 km (length) × 10 km (height). The existence of the near-surface conductors is consistent with the features of the observed data where the Φ 2 is higher than 45° at 20 Hz and low apparent resistivity data <100 Ωm is observed at a frequency of 100 Hz (Figure 2 and Figure S1-S6 in Supporting Information S1). The existence of deep conductor C1 is supported by the feature of the observed data where the Φ 2 is higher than 45° at 0.3 Hz around the northwest of the Takamatsu volcano ( Figure 2). Conductive zones extend from the east and west ends of deep conductor C1 to the near-surface conductors (Figure 3c). The conductive zone extending from the west end of C1 extends toward the Wasabizawa power plant (Figure 3c). Although a deep conductor is imaged at y = 9 km outside the coverage area of the MT measurements (Figure 3c), the MT data do not sufficiently constrain resistivity structures outside the coverage area of the measurements. Therefore, we exclude this conductor from our discussions in this study.  (Okada et al., 2014) are superimposed on the sections in (c and d). Note that depth is indicated as the depth below sea level, and x-and y-directions denote the geographic north and east, respectively. Black triangles denote MT sites within ±1 km of each profile line. Black dots denote the relocated hypocenters of earthquakes within ±150 m of each profile line (Okada et al., 2014). White lines denote the depth of a minimum of silica solubility in pure H 2 O under hydrostatic conditions (D Si_h : the depth is 1.3 km bsl and 1.9 km bgs at the drilling site N63-MS-6), as presented in Figure 6a.

of 21
The deep conductor C1 may be associated with a supercritical geothermal reservoir and melt. Therefore C1 is an important feature in the presented model. It is known that the inverse solution of MT data is non-unique, meaning that many models can fit the observed data equally well (Constable et al., 1987;Ishizu et al., 2021). We calculate the MT responses for models where the resistivity values of C1 are replaced by values between 1 and 100 Ωm (1, 3, 10, 30, and 100 Ωm) to determine how well the resistivity values are constrained by the data ( Figure S7 in Supporting Information S1). The RMS misfit is minimum for 10 Ωm and increases away from 10 Ωm ( Figure  S7 in Supporting Information S1). We measure the reliability of the resistivity values of C1 based on F values with a 95% confidence interval (Gresse et al., 2021;Yamaya et al., 2017). The threshold RMS misfit with a 95% confidence interval is 2.392. Forward models associated with C1: 6.5 Ωm < C1 < 35 Ωm obtain an RMS misfit <2.392. Forward models associated with low and high resistivities (i.e., <6.5 Ωm and >35 Ωm) show RMS values higher than F RMS = 2.392. In other words, these forward models are statistically rejected with a 95% confidence interval. Hence, we conclude that the resistivity range of C1 is 6.5 Ωm < C1 < 35 Ωm.

Shallow Hydrothermal System
First, a shallow geothermal system is investigated using the obtained resistivity model. The resistivity model reveals the near-surface conductors at depths between −1.3 and 0.5 km bsl (Figure 3). Similar near-surface conductors found in various volcanic areas have been interpreted as smectite-rich zones (Cherkose & Saibi, 2021;Kanda et al., 2019;Ledo et al., 2021;Pellerin et al., 1992;Tseng et al., 2020;Yoshimura et al., 2018). Drilling confirmed the smectite-rich zones in this field (New Energy and Industrial Technology Development Organization, 1990). Moreover, the distribution of near-surface conductors is consistent with the distribution of smectite-rich hydrothermal alteration zones at the surface (Nunohara et al., 2021). These results indicate that the near-surface conductors are smectite-rich zones. Smectite under moisture-rich conditions can be converted to illite at temperatures exceeding 200°C (Pytte & Reynolds, 1989;Wersin et al., 2007;Yamaya et al., 2013). Rocks containing illite have resistivity values 6-10 times higher than those of rocks with a similar proportion of smectite under the same temperature, porosity, and salinity conditions (Ussher et al., 2000). We set the boundaries of the smectite and illite zones at 50 Ωm because the depth at which the smectite-rich near-surface conductors reach a resistivity value of 50 Ωm is roughly consistent with the isotherm at 200°C (Figure 3). Hence, we interpret the resistive zone of >50 Ωm below the near-surface conductors as an illite-rich zone transformed from smectite above 200°C.
Conductive zones are found to extend from the east and west ends of deep conductor C1 to the near-surface conductors (Figure 3c). The conductive zone extending from the west end of C1 extends toward the Wasabizawa power plant (Figure 3c). The 3 He/ 4 He and 4 He/ 20 Ne ratios of the fumarolic gas at the Kawarage hot spring, located 3 km west of the Wasabizawa power plant, were 9.5 × 10 −6 and 180, respectively (Kita et al., 1992). This indicates that the helium at the Kawarage hot spring originates from magmatic gas. These results suggest that these conductive zones act as a path for magmatic fluid ascending from C1, resulting in anomalies in the helium ratios at the Kawarage hot spring. A conductor extending from the deep parts (25 km depth) to the surface at Naruko volcano also represents a path for ascending magmatic fluid (Ogawa et al., 2014).

Supercritical Geothermal Reservoir
The low resistivity of C1 suggests that it includes magmatic fluids and melt. C1 is imaged in the low Vp and Vs region, supporting that C1 contains magmatic fluids and melt ( Figure S8 in Supporting Information S1). The low Vp and Vs zones extend west compared to C1 ( Figure S8 in Supporting Information S1). We consider that the geometry difference largely comes from the model grid difference between the presented resistivity model and the velocity model (horizontal cell size is 300 m in the resistivity model and 6 km in the velocity model). The temperature of C1 is estimated to be above 400°C based on the temperature profile (Figure 4a). While few seismic events occurred inside C1, most occurred around the edge of C1 (Figure 3), indicating that seismic events may have resulted from fluid migration around the edge of C1. The depth of the seismicity cut-off is consistent with the top of C1, and it is known that the depth of the seismicity cut-off roughly corresponds to 400°C (Mitsuhata et al., 2001;Ogawa et al., 2001;Okada et al., 2014). The correspondence between the top of C1 and the seismicity cut-off also supports that C1 is >400°C.  en, 1958), suggests that the part of C1 at a temperature lower than 725°C is mostly crystallized and high conductivity values of the part of C1 <725°C is derived from the existence of magmatic fluids, rather than melt. A depth of 3.5 km bsl (or 4.1 km bgs at the drilling site N63-MS-6) corresponds to a temperature of 725°C and a pressure of 107 MPa for the lithostatic conditions ( Figure 4a). Therefore, we consider the upper part of C1 at a depth shallower than 3.5 km bsl as a potential supercritical geothermal reservoir (>400°C). The Vp/Vs ratio (Okada et al., 2014) is lower than 1.65 in the upper parts of C1 (Figure 3 and Figure S8 in Supporting Information S1). The low Vp/Vs ratio in the upper crust in NE Japan cannot be explained by the presence of melt but can be explained by the presence of a few volume percent of H 2 O with a large aspect ratio due to a higher bulk modulus of melt than a bulk modulus of H 2 O (Nakajima et al., 2001). Thus, the low Vp/Vs ratio in the upper parts of C1 also supports that the upper parts of C1 might contain a certain fraction of magmatic fluids. An increased Vp/Vs ratio is observed in the lower parts of C1, supporting the existence of melt in the lower parts of C1.
The typical NaCl content of magmatic fluids was estimated to be 5-10 wt%, based on fluid inclusion analysis of porphyry ore samples collected in various regions of the world (Heinrich, 2005). We assume that the magmatic fluid in our study area (subduction zone of NE Japan) also has the typical NaCl content of 5-10 wt%. The validity of the assigned 5-10 wt% NaCl content of magmatic fluids in our study area is supported by the result of the fluid inclusion analysis of Kofu Granite samples in central Japan (Kurosawa et al., 2010). Recent mantle xenolith studies reported the NaCl content of fluids in the mantle wedge beneath the Pinatubo volcano, Philippines (subduction zone) to be 5.1 wt% (Kawamoto et al., 2013). The reported value of 5.1 wt% NaCl content also supports that the assigned 5-10 wt% NaCl content is reasonable as an NaCl content of magmatic fluids in subduction zones. Phase states of NaCl-H 2 O fluids are different depending on the pressure, temperature, and salinity conditions (Afanasyev et al., 2018;Driesner & Heinrich, 2007;Sillitoe, 2010;Weis et al., 2012). We consider the phase states of the potential supercritical reservoir at a depth of 1.8 km bsl (temperature of 450°C) because shallow parts are easier to exploit than deeper parts. The pressure at the potential supercritical reservoir is  Sealing expected to be hydrostatic or lithostatic, or somewhere between the two (Figure 4b). We calculate the hydrostatic pressure assuming a water density of 1,000 kg/m 3 and lithostatic pressure assuming a rock density of 2,700 kg/ m 3 (Saishu et al., 2014).
The phase state of the supercritical fluid at 450°C and with a 5 wt% NaCl equivalent is a vapor-halite coexistence if the pressure is in the hydrostatic condition ( Figure 5). In the phase state of vapor-halite coexistence, solid halite and vapor coexist (Driesner & Heinrich, 2007). The supercritical fluid in the pore space exists as vapor with solid halite in this phase state. Although only a very small fraction of interstitial fluid in the halite can reduce electrical resistivity by orders of magnitude (T. Watanabe & Peach, 2002), the halite exists with vapor in this case. Therefore, we consider that the halite is dry. The resistivity of dry halite is reported as >10 8 Ωm (T. Watanabe & Peach, 2002). We assume that the vapor is resistive at 1,000 Ωm (Peacock et al., 2020). As a result, the resistivity of the vapor-halite coexistence in the pore space is estimated to be more resistive than 1,000 Ωm. Dry granite resistivity at a temperature of 450°C is 10 5 Ωm (Kariya & Shankland, 1983). These resistivity values suggest that the bulk resistivity of the supercritical geothermal reservoir with the phase state of vapor-halite coexistence is expected to be much more than 10 Ωm of C1. Thus, the phase state of vapor-halite coexistence would not apply to the supercritical geothermal reservoir.
The phase state of the supercritical fluid is liquid-vapor coexistence at 450°C if the pressure is below the average between hydrostatic and lithostatic conditions (e.g., 35 MPa; Figure 5). In the phase state of the liquid-vapor coexistence, the fluid exists as low-salinity vapor in equilibrium with a small fraction of hypersaline brine (Driesner & Heinrich, 2007). Laboratory measurement data reported that the conductivity of volcanic rock samples with steam in the pore space decreased by a factor of 20 compared to the conductivity of the rock samples with only liquid in the pore space (Milsch et al., 2010). In fact, the bulk resistivity of 100-1,000 Ωm has been reported for rocks with supercritical fluids of liquid-vapor coexistence (Gresse et al., 2021;Samrock et al., 2018), which are much more resistive than 10 Ωm of C1, indicating that the phase state of liquid-vapor coexistence would not also apply to the supercritical geothermal reservoir.  (Bakker, 2018). The magenta circle denotes the critical point (422°C and 34 MPa). Red, blue, and black crosses denote the pressure-temperature points at a temperature of 450°C (a depth of 1.8 km bsl) under lithostatic pressure, hydrostatic pressure, and average of lithostatic and hydrostatic pressure, respectively. (b) Isothermal pressure-composition sections for 450°C. The red and blue lines denote pressure at a depth of 1.8 km bsl under hydrostatic and lithostatic conditions. In the phase state of vapor-halite coexistence, solid halite and vapor coexist (Driesner & Heinrich, 2007). In the phase state of the liquid-vapor coexistence, the fluid exists as low-salinity vapor in equilibrium with a small fraction of hypersaline brine. In the F domain at pressures above the liquid-vapor surface, fluid properties can continuously vary between vapor-and liquid-like without any heterogeneous phase change (Driesner & Heinrich, 2007;Scott et al., 2015). The supercritical fluid in this domain is often called "single-phase fluid" (Driesner & Heinrich, 2007). We also use the term of single-phase fluid for the supercritical fluid in this region.
If the pressure is >42 MPa for the temperature of 450°C and 5 wt% NaCl content, no clear distinction between liquid and vapor is possible for the supercritical fluid in this condition (Driesner & Heinrich, 2007). Supercritical fluids in this state are often referred to as "single-phase fluids" (Driesner & Heinrich, 2007). In this article as well, supercritical fluids in this region will be referred to as "single-phase fluids". If the pressure is in the lithostatic condition, fluids in the supercritical geothermal reservoir exist as a "single-phase fluid" (Figure 5). A single-phase supercritical fluid >400°C is more electrically resistive compared to fluids at 300°C (Bannard, 1975;Kummerow et al., 2020;Nono et al., 2020). However, the conductivity of a single-phase supercritical fluid under the lithostatic pressure (64 MPa) and for 450°C and 5 wt% NaCl content shows a high value of 17.6 S/m even under a supercritical condition (N. Watanabe et al., 2021). This high conductivity (17.6 S/m) of a single-phase supercritical fluid under the lithostatic pressure implies that the supercritical fluid can reasonably explain 10 Ωm of C1. Therefore, we consider the fluid in the supercritical geothermal reservoir as a single-phase fluid under pressure close to the lithostatic pressure.
If a connection exists to the surface from the supercritical geothermal reservoir, the pressure should be under hydrostatic conditions in the supercritical geothermal reservoir and supercritical fluids move to the surface without the accumulation. Moreover, the above discussion suggests that pressures close to the lithostatic pressure are expected in the supercritical geothermal reservoir to explain the low resistivity (10 Ωm) of C1. Hence, sealing to weaken the connections should be considered as an explanation for the accumulation of supercritical fluids and the lithostatic pressure of the reservoir (Scholz, 2019;Sibson, 2020). Silica sealing has been found above potential supercritical reservoirs and may separate the hydrostatic and lithostatic regions (Fournier, 1999;Ingebritsen & Manning, 2010;Lowell et al., 1993;Manning, 1994;Saishu et al., 2014;Weatherley & Henley, 2013). To consider the potential of silica sealing in this field, we calculate the silica solubility using Loner AP software (Akinfiev & Diamond, 2009).
Two general mechanisms have been proposed for the occurrence of silica sealing above supercritical geothermal reservoirs. The development of silica sealing by the first mechanism occurs when upwelling magmatic fluids under lithostatic pressure encounter fluids under hydrostatic pressure (Fournier, 1999;Saishu et al., 2014). The first mechanism can be explained by the fact that at depths >1 km bsl, the solubility of silica in H 2 O under the hydrostatic pressure (black line in Figure 6a) is much smaller than the solubility under the lithostatic pressure (gray line in Figure 6a). When silica-rich magmatic fluids under lithostatic pressure encounter fluids under hydrostatic pressure, silica precipitates due to the difference in solubility. The development of silica sealing by the second mechanism occurs when meteoric water under the hydrostatic pressure moves downward (surface to deep part) driven by hydrothermal convection (Saishu et al., 2014). Silica solubility in water under the hydrostatic pressure increases with increasing depth, up to a depth of 1.2 km bsl (black line in Figure 6a). Then, the solubility decreases and has a minimum at a depth of D Si_h = 1.3 km bsl for the hydrostatic pressure indicating that silica precipitation can occur at a depth of D Si_h = 1.3 km bsl from downward-moving fluids (black line and green star in Figure 6a). D Si_h = 1.3 km bsl occurs at a temperature of 364.14°C and pressure of 18.44 MPa and coincides with the point of change from the liquid phase to the vapor phase for the downward-moving fluids (Figure 7). The solubility of silica in water depends on the density of the water (Tsuchiya & Hirano, 2007). The fluid flow from the higher density liquid region into the lower density vapor region causes a drop in silica solubility at D Si_h = 1.3 km bsl (Tsuchiya & Hirano, 2007).
To explain the silica precipitation caused by the first mechanism, the hydrostatic-lithostatic pressure boundary should exist before the first mechanism occurs. Therefore, we consider that the second mechanism occurred first, and the second mechanism would have caused the quartz precipitation at a depth of D Si_h = 1.3 km bsl and separates the hydrostatic and lithostatic regions at D Si_h . This interpretation is supported by correspondence between D Si_h and the depth of the top of C1 at the drilling point N63-MS-6. After the development of the silica sealing at D Si_h by the second mechanism, the first mechanism also contributes to enhancing the silica sealing at D Si_h . The silica-rich magmatic fluids under lithostatic pressure encounter fluids under the hydrostatic pressure at D Si_h , and the difference in silica solubility causes silica precipitation at D Si_h . Therefore, two mechanisms are currently contributing to the development of silica sealing at a depth of D Si_h = 1.3 km bsl (Figure 4c). The silica sealing also plays a role in trapping the upwelling supercritical fluids due to its low permeability.
We estimate the fluid fraction of a potential supercritical reservoir using the obtained resistivity model. The resistivity of dry granite is reported to be about 10 5 Ωm at a temperature of 400°C-550°C (Kariya & Shankland, 1983). The solid-phase resistivity for the supercritical geothermal reservoir is fixed at 10 5 Ωm. Reservoir pressure is assumed to be under the lithostatic condition ( Figure 4c). The resistivity of supercritical fluid at a depth of 1.8 km bsl (450°C and 64 MPa) are 17.6 and 36.0 S/m for 5 wt% and 10 wt% NaCl equivalent, respectively (N. Watanabe et al., 2021). The HS+ model suggests that 10 Ωm resistivity of C1 requires a fluid fraction of 0.9% and 0.4% for 5 wt% and 10 wt% NaCl equivalent, respectively ( Figure 8 and Table 1). Considering the resistivity range of C1 (6.5 Ωm < C1 < 35 Ωm), the fluid fraction of the potential supercritical reservoir at a depth of 1.8 km bsl is estimated to be 0.1%-1.3% (Figure 8 and Table 1). We also estimate the fluid fraction of the supercritical reservoir at a depth of 2.4 km bsl (550°C and 79 MPa). The resistivity of supercritical fluid is 5.4 and 11.5 S/m for 5 wt% and 10 wt% NaCl equivalent, respectively (N. Watanabe et al., 2021). Considering the resistivity range of C1 (6.5 Ωm < C1 < 35 Ωm), the fluid fraction of a potential supercritical reservoir at a depth of 2.4 km bsl is 0.4%-4.2% (Figure 8 and Table 1). Therefore, the fluid fraction of the supercritical reservoir is estimated to be 0.1%-4.2% with a salinity of 5-10 wt%. The 0.1%-4.2% fluid fraction of the supercritical reservoir estimated by the resistivity model is consistent with the porosity (1.7%) of the core sample of the Kakkonda granite near the bottom of the WD-1a well (Muraoka et al., 1998). A recent investigation of plagioclase alteration in mafic schists in NE Japan reported that significant nano-to micro-scale pores were generated in plagioclase during the alteration of 400°C-570°C, which resulted in the whole-rock porosity of 1.3% ± 0.2% (Nurdiana et al., 2021). The porosity values obtained by investigation of plagioclase alteration also support that our estimation of the fluid fraction of the supercritical reservoir is reasonable.

Notes.
We consider two different depths of the supercritical geothermal reservoir. The conductivity of supercritical fluids is calculated using a viscosity-dependent empirical model by N. Watanabe et al. (2021). Bulk resistivity of the supercritical geothermal reservoir represents the resistivity range of C1 determined by F values with a 95% confidence interval. The fluid fraction represents an estimated fluid fraction of the supercritical geothermal reservoir using the HS+ model. The resistivity of the solid phase is fixed at 10 5 Ωm (Kariya & Shankland, 1983).

Table 1 Summary of Fluid Fraction Estimation of the Supercritical Geothermal Reservoir
A petrologic study found dacitic and andesitic rocks in the Yuzawa geothermal field (Ban et al., 2007). Note that basaltic magmas can transform to andesite, dacite, and rhyolite through fractional crystallization (Wilson, 1989).
The dacitic melt has a shallower origin than the andesitic melt (Ban et al., 2007;Tatsumi et al., 2008). We assume that conductor C1 at a depth of 5 km bsl consists of dacitic melt and solid rock, and use the resistivity value reported by Laumonier et al. (2015) for dacitic melt. Table 2 describes the resistivity values of dacitic melt. Bulk resistivity is calculated using the Hashin-Shtrikman upper bound model (HS+; Hashin & Shtrikman, 1962), which predicts well-connected melt fractions (Naif et al., 2013). The resistivity of dry granite is reported to be about 1,000 Ωm at a temperature of 900°C-1,050°C (Kariya & Shankland, 1983). Thus, the solid-phase resistivity is fixed at 1,000 Ωm. The bulk resistivity of the dacitic melt and solid rock at 900°C and 150 MPa (at a depth of 5 km bsl) suggests that an 80% melt fraction with 2 wt% H 2 O is required to explain the 10 Ωm of C1 (Figure 9a). If we consider the lower bound of C1 (6.5 Ωm), the 100% melt fraction is considered for 2 wt% H 2 O (Figure 9a). However, these melt fractions are not reasonable. The water content in the dacitic melt can range up to 5.5 wt%, which is the water saturation condition (Wallace, 2005 Laumonier et al. (2015) for dacitic melt and those reported by Guo et al. (2017) for andesitic melt.  Figure 9. (a) The bulk resistivity of the dacitic melt and solid rock at 900°C and 150 MPa (a depth of 5 km bsl) as a function of melt fraction and water content in the dacitic melt. The resistivity of the dacitic melt is calculated based on Laumonier et al. (2015). Solid rock resistivity is fixed at 1,000 Ωm (Kariya & Shankland, 1983). Bulk resistivity is calculated using Hashin-Shtrikman upper bound model (HS+). (b) The bulk resistivity of the andesitic melt and solid rock at 1,050°C and 250 MPa (a depth of 9 km bsl) as a function of melt fraction and water content in the andesitic melt. The resistivity of the andesitic melt is calculated based on Guo et al. (2017). The gray zones denote the resistivity range with a 95% confidence interval for C1 based on the F-test (6.5 Ωm < C1 < 35 Ωm). a 20% melt fraction can explain 10 Ωm (Figure 9a). Even if we consider the lower bound of C1 (6.5 Ωm), a 30% melt fraction is considered for 5 wt% H 2 O. Therefore, we consider that water content in the dacitic melt is 5 wt% H 2 O. The upper bound of C1 (35 Ωm) requires a 6% melt fraction for 5 wt% H 2 O. These findings suggest that C1 at a depth of 5 km bsl is a dacitic melt with 5 wt% H 2 O and a melt fraction of 6%-30%.
The andesitic melt is expected to be located in the deep zones of this field (Ban et al., 2007). The conductor C1 at a depth of 9 km bsl is assumed to consist of andesitic melt and solid rock. The resistivity of the andesitic melt is calculated according to Guo et al. (2017). Table 2 describes the resistivity values of andesitic melt. The water content in the andesitic melt can range up to 5.5 wt% for 1,100°C-1,300°C and 200 MPa (Botcharnikov et al., 2006). Water is enriched in the melt during fractional crystallization, and therefore the dacitic melt should have a higher water content than the andesitic melt. This suggests that the water content in the andesitic melt is lower than the estimated value of 5 wt% H 2 O in the dacitic melt. To explain 10 Ωm of C1 at 1,050°C and 250 MPa (at a depth of 9 km bsl), the melt fraction is estimated to be 20% and 8% for 2 wt% and 4 wt% H 2 O contents, respectively ( Figure 9b). The lower bound of C1 (6.5 Ωm) requires 30% and 12% melt fractions for 2 wt% and 4 wt% H 2 O contents, respectively. The upper bound of C1 (35 Ωm) requires 6% and 2% melt fractions for 2 wt% and 4 wt% H 2 O content, respectively. Hence, our resistivity model suggests that C1 at a depth of 9 km bsl consists of andesitic melt with a melt fraction of 2%-30% and 2-4 wt% H 2 O in the melt.

Mechanism of Evolution of a Supercritical Geothermal Reservoir
Our MT data reveal a potential supercritical geothermal reservoir, melt, and a shallow geothermal system in the Yuzawa geothermal field in NE Japan (Figure 3). Based on our resistivity model, we propose a mechanism for the evolution of a supercritical geothermal reservoir ( Figure 10). The dacitic and andesitic melt below the supercritical geothermal reservoir may supply magmatic fluids to the supercritical geothermal reservoir (Blundy et al., 2021;Heinrich, 2005;Sillitoe, 2010). Upwelling supercritical fluids supplied from the melt are found to become trapped under a less-permeable silica sealing, and supercritical fluids accumulated below the silica sealing as a result ( Figure 10). Silica sealing separates the hydrostatic and lithostatic regions (Figure 4c). Hence, the supercritical geothermal reservoir below the silica sealing is under lithostatic pressure, and fluids in the supercritical geothermal reservoir are interpreted as single-phase supercritical fluids. The supercritical geothermal reservoir is imaged as a low resistivity anomaly of C1 ( Figure 3). The fluid fraction of the supercritical reservoir is estimated to be 0.1%-4.2% (Figure 8). Episodic supplies of magmatic fluids from the melt increase the pressure of the supercritical geothermal reservoir, with high levels of pressure subsequently breaking the silica sealing. Fluids leaking from the silica sealing moved to the surface ( Figure 10). The fluids are detected as conductive zones extending from the deep conductor C1 to the near-surface conductors (Figure 3c). The smectite-rich zones near the surface act as a cap layer owing to their low permeability (Revil et al., 2019). The upwelling fluids that leak from the silica sealing mixes with the meteoric water during the ascent, and the mixed hydrothermal fluids are trapped below the cap layer. Power plants in operation in the area use hydrothermal fluids (<300°C) trapped below the impermeable cap of smectite to generate electric energy (New Energy and Industrial Technology Development Organization, 1990). Since our MT data reveal a potential supercritical geothermal reservoir in this geothermal field, power plants in this field potentially may use supercritical geothermal reservoirs to increase their power generation in the future.

Conclusions
Supercritical geothermal reservoirs are next-generation energy resources that yield higher productivity than those obtained from conventional geothermal reservoirs at temperatures <350°C. Although understanding the fluid fraction and spatial distribution of supercritical geothermal reservoirs Figure 10. (a) West-east cross-section at x = 0.5 km of the inverted resistivity model. This cross-section is the same as Figure 3c. Symbols are the same as in Figure 3. (b) Schematic model of a geothermal system inferred from our resistivity model.
is necessary for their assessment as next-generation geothermal resources, these characteristics remain poorly understood. In this study, we applied the MT method to the Yuzawa geothermal field in NE Japan to estimate the fluid fraction and spatial distribution of a potential supercritical geothermal reservoir. Our main findings can be summarized as follows: 1. The MT data reveal a supercritical geothermal reservoir, melt, and shallow geothermal system in the Yuzawa geothermal field 2. The supercritical geothermal reservoir (>400°C) with a size of 3 km (width) × 5 km (length) is located at a depth of 2.5-6.0 km below the surface. The fluid fraction of the supercritical reservoir is estimated to be 0.1%-4.2% with salinity of 5-10 wt% 3. Silica sealing may exist above the potential supercritical geothermal reservoir, separating the hydrostatic and lithostatic regions. The supercritical geothermal reservoir is under lithostatic pressure, and fluids in the reservoir exist as single-phase supercritical fluids 4. Dacitic and andesitic melt exist below the supercritical geothermal reservoir. The melt supplies magmatic fluid to the supercritical geothermal reservoir 5. We develop a mechanism to explain the evolution of a supercritical geothermal reservoir, wherein upwelling supercritical fluids supplied by the melt are trapped under less permeable silica sealing and accumulated there 6. The supercritical fluid breaking from the silica sealing provides upward-moving fluids to the surface

Data Availability Statement
The MT data used in this study (Ogawa et al., 2011(Ogawa et al., -2020 were archived at the Data Management Center of the Incorporated Research Institutions for Seismology SPUD EMTF database (http://ds.iris.edu/spud/emtf). The following two formats are available for all MT data: electrical data interchange (EDI) format and electromagnetic transfer function extensible markup language (EMTF XML) format. The EDI format is a standard MT data format constructed by the Society of Exploration Geophysicists (Wight, 1988) and the description of EDI format is provided at the following website: https://www.seg.org/Portals/0/SEG/News%20and%20Resources/Techni-cal%20Standards/seg_mt_emap_1987.pdf. The EMTF XML format is a new, self-describing, searchable, and extensible way to store MT data (Kelbert, 2020). The open-source Python package MTpy can be used to analyze the MT data (Kirkby et al., 2019;Krieger & Peacock, 2014). The quartz solubility, pressure, and temperature data for Figure 6 are available on the Zenodo open-access repository operated by the CERN: https://doi.org/10.5281/ zenodo.5952138. One of the authors (K. Ishizu) would like to acknowledge the support from the Hiki Foundation, Tokyo Institute of Technology. The authors thank W. Siripunvaraporn for allowing us to use WSINV3DMT, and T. Okada for providing seismic data and velocity models. The authors are grateful to N. Watanabe for helpful discussions on the resistivity of supercritical fluids, Y. Kawada for meaningful discussions on supercritical fluids, R. Yamada for his help with the geological interpretation of the resistivity model and A. Kelbert for her help in archiving our MT data. The authors would like to thank two anonymous reviewers, the associate editor M. Moorkamp and the editor D. Schmitt for their detailed and helpful comments that led to significant improvements in this article.