grace data recorded postseismic increase of the gravity field by 6 galviscoelastic relaxation and afterslip were examined with the gravity datathe grace data constrain the biviscous rheology model for the asthenosphere grace data recorded postseismic increase of the gravity field by 6 gal viscoelastic relaxation and afterslip were examined with the gravity data the grace data constrain the biviscous rheology model for the asthenosphere
the great ( mw 9.1 ) 11 march 2011 tohoku - oki earthquake disturbed not only the ground and ocean surface but also the earth 's gravity field permanently .
it generated transient seismic and tsunami waves and produced large coseismic stress change that led to gradual stress redistribution and subsequent postearthquake crustal motions .
the dense geodetic data such as gps , interferometric synthetic aperture radar , and acoustic seafloor measurements reported coseismic displacement ( a few tens of meters of horizontal motion and a few meters of vertical motion ) and rapid postseismic displacement that yielded 10% of the equivalent seismic moment of the main shock in 2 weeks [ ozawa et al . , 2011 ; sato et al . , 2011
additionally , the gravitational effect of the earthquake has been measured through changes of instantaneous relative motions of two gravity recovery and climate experiment ( grace ) satellites coorbiting around 500  km altitude [ han et al . ,
subsequently , the time series of monthly grace gravity field data have been processed by optimizing the signal over the earthquake region [ matsuo and heki , 2011 ; cambiotti and sabadini , 2012 , 2013 ; wang et al .
most of the large earthquakes are followed by seismic or aseismic sliding ( afterslip ) on the coseismic fault and its lateral extension , by viscoelastic flow as a process of gradual relaxation of coseismic stresses , and/or by deformation induced by pore fluid migration in response to coseismic stresses [ e.g. , cohen , 1999 ; wang et al .
inverted the postseismic gps displacement data over the 1.5  years after the 2011 tohoku - oki earthquake for afterslip and viscoelastic relaxation models , and concluded afterslip is a dominant process while viscoelastic relaxation is responsible for only 11% of the postseismic moment .
however , they did not evaluate biviscous processes that may explain transient deformation during the first few months as well as afterslip .
for example , after the 2004 sumatra - andaman earthquake , in addition to gps displacement data at teleseismic distances [ pollitz et al . , 2006 ] , the grace gravity data were useful for constraining postseismic processes consistent with biviscous relaxation [ han et al . , 2008 ; panet et al . , 2010 ; hoechner et al . , 2011 ] .
it is imperative to constrain the complex response of the earth to the 2011 tohoku - oki earthquake to understand the rheology of the crust and mantle and the cycle of strain accumulation and release , ultimately , in order to estimate the subsequent seismic hazards in the surrounding region [ wang et al . , 2012a ] . to this end
, we analyzed spatially and temporally continuous grace gravity observations implying the broadscale ( 500  km ) deformation of the surface and interior in response to regional - scale postearthquake stress / strain redistribution and evaluated two alternate postseismic processes of viscoelastic relaxation and afterslip .
additionally , we calculated the high - resolution gravimetric response to this event and used these results with other available geodetic data to further discriminate between alternate mechanisms .
we analyzed a total of 11  years of grace gravity data from 2003 to 2013 using the rl05 level-2 ( l2 ) monthly data products generated by center for space research , university of texas [ tapley et al .
, 2004 ] . the march 2011 data were not used since the earthquake occurred in the middle of that month .
the spherical harmonic coefficients up to degree and order 40 were used to represent gravity variations at a spatial resolution of 500  km , which is a conservative estimate of the grace data resolution .
we first determined an analytical model of annual and semiannual sinusoids and a linear trend in the grace data by fitting the model to the monthly grace time series from january 2003 to february 2011 ( i.e. , prior to the earthquake ) .
we extended the estimated trend out to the end of 2013 and removed it from the entire set of grace data before the analysis .
figure  1a presents the grace gravity changes calculated by subtracting the mean gravity field over a preearthquake period ( february 2008 to february 2011 ) from the mean field over a postearthquake period ( april 2011 to december 2013 ) .
therefore , it illustrates the coseismic deformation as well as the postseismic deformation averaged over a period from 11 march 2011 to december 2013 .
next , we computed the average coseismic gravity changes ( figure  1b ) using the seismic centroid moment tensor ( cmt ) solutions ( global cmt and u.s .
geological survey w - phase moment tensor ) and two seismic finite fault models , models ii and iii , from shao et al . .
the finite fault models we used present the largest slip updip of the hypocenter close to the trench consistent with the recent analysis including tsunami and seafloor data presenting the rupture out to the trench [ lay et al . , 2011 ; iinuma et al . ,
the dominant negative anomaly found in the seismic models indicates extensive crustal density decrease by volume extension under the influence of the ocean mass redistribution [ han et al . , 2006 , 2013 ;
( a ) the grace observed gravity change following the 2011 tohoku - oki earthquake computed by differencing the 3 year mean field before and after the earthquake .
( b ) the average synthetic gravity change caused by the coseismic deformation computed using four seismic models .
( c ) the postseismic grace gravity change computed by subtracting the seismic model prediction shown in figure  1b from the mean difference field shown in figure  1a .
( d ) the same as figure  1a but during the 3 years between 20082010 and 20052007 ( no large earthquake in - between ) .
figure  1c elucidates the postseismic gravity change computed by removing the coseismic ( elastic ) contribution ( figure  1b ) from the grace data ( figure  1a ) .
it indicates that the magnitude of postseismic deformation is a fraction of the coseismic deformation . during the first 3  years after the main shock , on average ,
6 gal of postseismic gravity change was estimated mostly around the epicenter . during a time period that does not include the perturbation due to a large earthquake ( for example
, the 20082010 mean field minus the 20052007 mean field ) , the mean field difference is within 2 gal ( figure  1d)a measure of grace data noise and inherent gravity variations in this region .
while the coseismic gravity change produced the positive anomaly offshore and the stronger negative anomaly landward , the postseismic gravity change yielded the positive anomaly prominently around the epicenter where the coseismic perturbation was small . to understand this
, we further examine the temporal dependence of the postseismic gravity data and compare it with models of viscoelastic relaxation and afterslip .
the temporal characteristics of postseismic gravity changes provide constraints on the physical processes and the rheological structure of the earth .
in particular , pollitz [ 1997b ] found that asthenospheric viscosity controls the temporal pattern of long - wavelength postseismic relaxation following a large megathrust event . using the global normal mode formalism and viscous relaxation theory given in pollitz [ 1997b ] and the fault geometry notation of kanamori and cipar [ 1974 , figure  10 ] and kanamori and given [ 1981 , figure  1 ]
, we express the gravitational potential change induced by earthquake as follows : 
 1where a point dislocation source is at the ( north ) pole and described with the moment tensor components , mrr , mr , mr , m    m , and , m. the five independent and orthogonal basis functions of the earthquake geopotential change , gmrr , gmr , gmr , gm    m , and gm , are excited by each moment tensor component with a physical dimension of m / s per nm .
they are computed with the eigenfunctions of viscoelastic normal modes ( pollitz [ 1997b ] , pollitz et al . , and han et al .  
they are dependent only on the earth 's viscoelastic structure ( i.e. , precomputable regardless of seismic sources ) .
the examples of the spatial pattern of each function for the elastic response ( evaluated at t  =  0 ) are shown in figure  2 ( left ) .
( left ) the spatial patterns of the gravity changes to be excited by five independent moment tensor components of mrr , mr , mr , m    m , and , m , respectively .
the exact scale is dependent on the moment tensor ( the negative to positive values are depicted with blue to red colors ) .
the coordinates of the dislocation source are 38.5n , 142.6e ( depicted as a red star ) , and its depth is 20  km .
( right ) the corresponding temporal variations of each component of the gravity change before and after the earthquake .
the mean value of the synthetic gravity changes computed from various seismic models is shown in thick black solid line , and the 1 sigma variation of the models is depicted in thin black dashed line . depicted
are viscoelastic gravity changes computed from the maxwell asthenosphere model with the viscosity of 5    10  pa  s ( blue ) , 10  pa  s ( green ) , 5    10  pa  s ( red ) , and 10  pa  s ( cyan ) and from the biviscous ( burgers body ) asthenosphere model with transient viscosity of 10  pa  s and steady state viscosity of 10  pa  s ( purple ) .
the temporal evolution of each gravity component , mrr , mr , mr , m    m , and m , are simply obtained by the inner product of each basis function ( e.g. , gmrr ) and the grace gravity data over unit sphere ( see the supporting information for the explicit forms ) .
the observed temporal changes from grace are presented in figure  2 ( right ) .
the gravity changes corresponding to the moment tensor components of mr and m can be resolved with the highest accuracy of 2030    10  nm ( or mw = 8.28.3 ) , while the other components have the accuracy of 80    10  nm ( or mw = 8.6 ) .
the north - south sampling nature of the grace intersatellite tracking enables the best sensitivity in estimation of mr and m components , 34 times better than other components .
this indicates the grace 's sensitivity to earthquakes as small as mw = 8.2 in this region , depending on the focal mechanism .
the average and standard deviation of four seismic models are shown as thick black solid and thin black dashed lines in figure  2 , respectively .
the first few months of grace data after the earthquake agree with the coseismic changes from seismic models ; however , subsequently , there are substantial variations observed particularly in mrr and m    m. there is no significant postseismic change discerned from other components .
we used a spherically stratified viscoelastic earth model that includes a global ocean layer ( 3  km thick ) , an elastic lithosphere ( 60  km thick ) overlying an asthenosphere extending to a depth of 220  km , and an upper ( 220670  km ) and lower ( 6702891  km ) mantle [ pollitz et al . , 2006 ] .
the elastic thickness of 63  km near the japan trench was estimated independently from the marine gravity and flexural analysis [ levitt and sandwell , 1995 ] .
the stratified density structure is consistent with preliminary reference earth model [ dziewonski and anderson , 1981 ] .
the upper and lower mantle was modeled with a maxwell rheology with the viscosity of 10  pa  s and 10  pa  s , respectively .
we tested five different rheological models for the asthenosphere ; ( 1 ) maxwell viscosity of 5    10  pa  s , ( 2 ) 10  pa  s , ( 3 ) 5    10  pa  s , ( 4 ) 10  pa  s , and finally , ( 5 ) transient ( kelvin ) viscosity of 10  pa  s and steady state ( maxwell ) viscosity of 10  pa  s with biviscous ( burgers body ) rheology [ ivins and sammis , 1996 ] .
the lithosphere overlying the asthenosphere is assumed to be elastic by assigning a very large viscosity of 10  pa  s. we use the computer code visco1d [ pollitz , 1997b ] and its modification to compute the postseismic gravity response to a dislocation source for alternate spherical viscoelastic models [ pollitz , 1997a ; pollitz et al . ,
the pollitz [ 1997b ] method is applicable to any deformation field represented as a viscoelastic normal mode sum on a laterally homogeneous spherical model , specifically , a weighted spherical harmonic sum of displacement - stress vectors .
although originally designed for the case of static gravity changes , it is also directly applicable to the case of postseismic gravity changes by employing equation ( 21 ) of pollitz [ 1997a ] as the viscoelastic normal mode sum .
figure  2 presents the various synthetic viscoelastic gravity changes , depicted in solid colored curves .
these viscoelastic simulations were computed with the cmt solution obtained from han et al . .
the equivalent simulations with a finite fault model such as model ii of shao et al .   are shown in the figure s1 in the supporting information , and they are not different significantly at the spatial scale of consideration because of the concentrated coseismic slip of this earthquake .
for the components mrr and m    m , the model runs with the asthenospheric maxwell viscosity of 5    10 and 10  pa  s ( blue and green curves , respectively ) predict rapid gravity change within a few months after the main shock .
this may be reconciled with the short - term grace observations ; however , the long - term changes with these viscosities are predicted to be overly large . the simulations with the asthenospheric maxwell viscosity of 5    10 and 10  pa  s ( red and cyan curves , respectively )
do not reproduce the rapid change observed during the first few months but capture the cumulative long - term change .
the rapid change followed by gradual variation in the grace gravity data are best reproduced with the biviscous relaxation model with the transient viscosity of 10  pa  s and steady state viscosity of 10  pa  s ( purple ) . for the other components of mr , mr , and m that do not show discernible postseismic gravity change observation , any model with the steady state viscosity smaller than 5 
overall , the biviscous asthenospheric model satisfactorily explains the postseismic gravity changes observed by grace , similar to the findings for the 2004 sumatra - andaman earthquake [ pollitz et al .
, 2006 ; han et al . , 2008 ; panet et al . , 2010 ; hoechner et al . , 2011 ]
alternatively , the afterslip on the coseismic rupture surface and its downdip extension is also a viable candidate to explain postseismic deformation after large megathrust earthquakes .
we use the finite fault models of afterslip to quantify the predicted gravity for comparison with the time series of grace gravity data .
the ( kinematic ) afterslip model was constrained by gps measurements of postseismic displacement accumulated over various intervals ( 7 , 20 , 48 , 107 , 200 , 289 , 381 , and 564  days ) up to the 1.5  years after the main shock , assuming that the postseismic displacement is entirely due to afterslip [ diao et al . ,
the maximum cumulative slip of 4  m over the period of 1.5  years and the equivalent seismic moment of 2.3    10  nm was found by elastic dislocation modeling of postseismic gps data [ ozawa et al . , 2011 ; diao et al . ,
, we calculated the deformation associated with afterslip using dislocation theory [ e.g. , cohen , 1999 ] . figure  3a presents the grace postseismic gravity data ( same as figure  1c ) , and three different locations where the positive and negative ( to a lesser extent ) gravity changes are identified . the gravity change of the mrr component is responsible mostly for the positive gravity at p2 , while the one of m    m is for the negative gravity change at p1 and p3 . the grace gravity time series at each location are presented with the error estimates of 2 gal in figure  3b . the elastic gravity change was computed with model ii of shao et al .
, and the viscoelastic gravity change was computed using the same seismic source model and the biviscous earth 's model characterized with the transient viscosity of 10  pa  s and steady state viscosity of 10  pa  s for the asthenosphere .
the time series of the afterslip gravity change from the finite fault model of diao et al
( a ) grace observation of postseismic gravity change ( same as figure  1c ) . ( b ) time series of the grace gravity change at p1 , p2 , and p3
. an arbitrary bias of 17 and 17 gal was added to the ones of p1 and p3 , respectively , for plotting clarity .
the model prediction of coseismic slip , afterslip , and viscoelastic relaxation is shown , respectively , as black , blue , and red lines .
( c ) a snapshot of postseismic gravity change 1 year after the main shock from the viscoelastic relaxation model .
( d ) same as figure  3c but from the afterslip model . the snapshot of the synthetic gravity change 1 year after the main shock is depicted in figures  3c and 3d , respectively , for the viscoelastic relaxation and the afterslip models .
both postseismic processes produce similar gravity change in spatial and temporal pattern as well as the magnitude .
the positive anomaly around the epicenter ( p2 ) is a dominant feature consistently found in the grace measurements and in both viscoelastic relaxation and afterslip models .
the secondary anomalies with 34 times smaller magnitude are the negative changes in the west and the east of the positive anomaly ( p3 and p1 , respectively ) .
the afterslip model predicts a slightly larger anomaly in the west than the one in the east .
so far , we examined the postseismic  free - air  gravity changes measured from grace and computed from viscoelastic relaxation and afterslip models at a spatial resolution of 500  km ; it is difficult to discriminate these competing postseismic processes from the grace data alone for the period of a couple of years after the main shock . in order to further understand gravity change and deformation associated with these postseismic mechanisms
, we computed the gravity effect of surface vertical deformation ( uplift / subsidence of the ocean - seafloor interface ) and interior deformation ( bouguer gravity anomaly ) separately [ han et al . , 2006 ] .
the gravity change of seafloor vertical motion by viscoelastic relaxation ( figure  4a ) is characterized by a dominant positive anomaly around the epicenter implying broadscale ( 500  km ) seafloor uplift .
the bouguer anomaly of viscoelastic relaxation ( figure  4b ) , computed by subtracting the effect of seafloor vertical motion from the free - air anomaly , shows a similar pattern of gravity change but with a magnitude 23 times smaller .
the bouguer gravity change is explained mostly by vertical deformation of the density interfaces of the upper and lower crust and of the crust and upper mantle ( moho ) .
it indicates that the broadscale uplift gradually occurs with negligible change in bulk volume ( or density ) as viscous asthenosphere and mantle yield to the coseismic stresses .
( a , c , and e ) gravity change due to seafloor vertical motion computed , respectively , from biviscous viscoelastic relaxation , afterslip , and coseismic slip .
( b , d , and f ) same as figures  4a , 4c , and 4e , respectively , but for the bouguer gravity . in comparison to the uplift by viscoelastic relaxation ,
the afterslip model predicts twice larger seafloor uplift around the epicenter as highlighted by the doubled magnitude of the positive gravity ( figure  4c )
. the bouguer gravity of afterslip ( figure  4d ) is characterized by a negative anomaly orthogonal to the strike direction .
this negative anomaly is primarily due to volume expansion or dilatation associated with the downdip afterslip within the lower crust and the upper mantle .
we also computed the gravity changes of seafloor vertical motion and bouguer anomaly ( figures  4e and 4f , respectively ) from the coseismic slip model , model ii of shao et al . .
the majority of the coseismic slip is estimated at the top 25  km of the crustal layer , while the afterslip is concentrated mostly at depths of 4050  km [ ozawa et al .
deeper sources like the downdip afterslip yield smaller volume change due to larger rigidity and bulk modulus at the mantle depth , while shallower coseismic slip associated with smaller rigidity and bulk modulus at the crustal depth yields larger volume change . as a result ,
the deep afterslip produces the smaller bouguer gravity anomaly ( figure  4d ) relative to the gravity from seafloor vertical deformation ( figure  4c ) and thus results in the smaller negative change in the free - air gravity ( figure  3d ) . on the other hand , the coseismic bouguer gravity anomaly
( figure  4f ) is as large as the coseismic vertical displacement ( figure  4e ) , which reduces the positive contribution to the free - air gravity ( figure  1b ) .
the downdip afterslip also moves the peak gravity change landward ( west ) , relative to the coseismic gravity change .
we repeated the same calculations of synthetic gravity changes for the viscoelastic , afterslip , and coseismic deformation but at a resolution of 55  km by extending the spherical harmonic degree and order up to 360 and by using a finite fault model shao et al .  
we found that the high degree expansion of the viscoelastic simulation predicts local subsidence superimposed on the broadscale uplift across the japan trench .
while the east coast of the island and the east of the trench undergo uplift , the area localized above the rupture zone undergoes subsidence , as implied by the positive and negative gravity change , respectively . this is the very characteristic pattern of surface vertical motion followed by thrust earthquakes when the elastic lithosphere overlying the viscoelastic asthenosphere is only  partially  ruptured [ melosh , 1983 ; cohen , 1999 ] .
melosh   found that the top boundary of the viscous asthenosphere experiences local extension beneath the fault if the fault stops within the lithosphere , which causes the local subsidence on top of regional uplift .
however , if the fault cuts through the lithosphere , the boundary of the asthenosphere is subject to net shortening and the surface undergoes uplift .
sato et al . , burgmann et al . , and hu et al .   reported the postseismic subsidence of the seafloor above the rupture zone by a few decimeters within a couple of years after the main shock using the gps - acoustic positioning and seafloor pressure data , while the gps data implied postseismic uplift along the coast by a few centimeters [ japan coast guard , 2012 ; japan coast guard and tohoku university , 2013 ] .
such subsidence may yield tens of gal gravity decrease ( 7 gal per 10  cm of seafloor subsidence ; the bouguer correction with the density contrast between the crust and ocean ) .
no seafloor measurement is available in the pacific plate , to the east of the trench , where viscoelastic uplift is predicted .
the local subsidence from the viscoelastic model is diminished by the uplift in the periphery of the overthrust and underthrust blocks at a broadscale resolution , which is what the grace data observed .
does not produce subsidence , although the model is consistent with onland postseismic gps motions .
reported such inconsistency in the postseismic vertical motion from afterslip following the great 2006 and 2007 kuril earthquakes .
hu et al .   found that the poroelastic rebound may yield uplift above the rupture zone by a few decimeters in 2 years .
although all of these postseismic processes are likely present and compensate each other at times , we hypothesize the viscoelastic relaxation triggered by the partially ruptured elastic lithosphere is a main driver of the local subsidence above the rupture region .
grace observed postseismic gravity increase by 6 gal ( 4050% of the coseismic gravity change ) within a couple of years after the 2011 tohoku - oki earthquake .
the postseismic free - air gravity changes from viscoelastic relaxation and afterslip are similar at the scale of 500  km commensurate with the grace data resolution . however ,
the bouguer gravity responses are distinct ; the bouguer anomaly of viscoelastic relaxation is positively correlated with the free - air anomaly , while the bouguer anomaly of afterslip is negatively correlated . the broadscale postseismic gravity observation and the local subsidence measurement above the rupture zone
the viscosity is the most important factor to govern evolution of the long - wavelength stress field after the earthquake .
large - scale monitoring of the postseismic deformation ( gravity as well as vertical displacement ) will help characterize the rheological properties of the earth 's interior that may lead to a better understanding of increased seismic hazard following great earthquakes .
 the analysis of grace gravity data revealed postseismic gravity increase by 6 gal over a 500 km scale within a couple of years after the 2011 tohoku - oki earthquake , which is nearly 4050% of the coseismic gravity change . 
 it originates mostly from changes in the isotropic component corresponding to the mrr moment tensor element . 
 the exponential decay with rapid change in a year and gradual change afterward is a characteristic temporal pattern . 
 both viscoelastic relaxation and afterslip models produce reasonable agreement with the grace free - air gravity observation , while their bouguer gravity patterns and seafloor vertical deformations are distinctly different . 
 the postseismic gravity variation is best modeled by the biviscous relaxation with a transient and steady state viscosity of 1018 and 1019  pa  s , respectively , for the asthenosphere . 
 our calculated higher - resolution viscoelastic relaxation model , underlying the partially ruptured elastic lithosphere , yields the localized postseismic subsidence above the hypocenter reported from the gps - acoustic seafloor surveying . 