Solar wind reconstruction from magnetosheath data using an adjoint approach

We present a new method to reconstruct solar wind conditions from spacecraft data taken during magnetosheath passages, which can be used to support, e.g., magnetospheric models. The unknown parameters of the solar wind are used as boundary conditions of an MHD (magnetohydrodynamics) magnetosheath model. The boundary conditions are varied until the spacecraft data matches the model predictions. The matching process is performed using a gradientbased minimization of the misfit between data and model. To achieve this time-consuming procedure, we introduce the adjoint of the magnetosheath model, which allows efficient calculation of the gradients. An automatic differentiation tool is used to generate the adjoint source code of the model. The reconstruction method is applied to THEMIS (Time History of Events and Macroscale Interactions during Substorms) data to calculate the solar wind conditions during spacecraft magnetosheath transitions. The results are compared to actual solar wind data. This allows validation of our reconstruction method and indicates the limitations of the MHD magnetosheath model used.


Introduction
The interaction of the solar wind with a planet and its magnetic field, such as the Earth's internal dipole field, causes the current system of the magnetosphere.Thereby, currents in the outer region of the magnetosphere, the magnetosheath, are related to the deflection of the solar wind plasma around the planet.These currents, far from the planetary surface, can be observed by a spacecraft crossing the magnetosheath with its boundaries being the bow shock and magnetopause.The influence of the solar wind on the currents in the inner magnetosphere can also be directly observed by ground magnetometer data.Modifications of the strength of these currents can be expressed by the indices of geomagnetic activity depending on the solar wind conditions.This dependence was used by Andreasen (1997) to reconstruct average solar wind conditions and their gross temporal variations in the pre-spacecraft era.Thereby, solar wind spacecraft observations were related to the indices of geomagnetic activity using a multi-regression method.Then, the solar wind conditions were estimated from ground magnetometer data back to 1926.A different approach to relate solar wind conditions to the indices of geomagnetic activity was introduced by Kondrashov et al. (2014).They use a singular spectrum analysis to fill gaps of the solar wind data in the period of 1972-2013.This analysis estimates hourly solar wind conditions during periods where no in situ measured solar wind information is available.It was shown that the reconstructed solar wind information can improve empirical models such as the Tsyganenko model (Tsyganenko and Sitnov, 2005).
In contrast to the previous approaches, Nagatsuma et al. (2015) included spacecraft observations to reconstruct solar wind conditions.During a geomagnetic storm in 1989, the magnetopause was near the geostationary orbit due to the extremely high solar wind dynamic pressure.The magnetopause location was determined by the Geostationary Operational Environmental Satellite (GOES).This magnetopause location was related to the solar wind conditions with the C. Nabert et al.: Solar Wind Reconstruction magnetopause model by Shue et al. (1998).This reconstruction approach using the spacecraft data with a high time resolution allows one to estimate the solar wind conditions more precisely on a much higher time resolution compared to the previous approaches.
We extend this approach by using not only the magnetopause location observed by a spacecraft, but by using spacecraft data obtained everywhere in the magnetosheath.Analogously to the previous studies, the reconstruction can be used to provide solar wind information if in situ solar wind observations are not available.This determines the solar wind parameters, which are the drivers of the magnetospheric current system and are usually required in magnetospheric and ionospheric models (e.g., Nagata et al., 2008;Korth et al., 2010;Zhang et al., 2013;Baker et al., 2013) We reconstruct the solar wind conditions during a spacecraft magnetosheath passage using the magnetosheath model introduced by Nabert et al. (2013).The solar wind parameters are boundary conditions of the magnetosheath model and they are varied until the model predictions match the magnetosheath data.The corresponding adjoint magnetosheath model is derived to perform this time-consuming procedure.Adjoint models are broadly used in fluid dynamic optimization problems such as drag minimization (e.g., Jameson, 1988;Othmer, 2008Othmer, , 2014;;Meader and Martins, 2012), or in seismology (e.g., Fichtner et al., 2006).We will introduce and discuss the problems of how to deal with a magnetosheath model when an adjoint approach is used.To verify our method, we use data from the THEMIS (Time History of Events and Macroscale Interactions during Substorms) mission (Angelopoulos, 2008) to reconstruct the solar wind conditions at Earth's orbit.The reconstructed conditions are compared to OMNI solar wind monitor data (see http://omniweb.gsfc.nasa.gov/)as well as THEMIS data from a second spacecraft in the solar wind.This indicates the limitations of the magnetosheath model used.

The forward magnetosheath model
For the solar wind reconstruction, we chose the zeroth order MHD (magnetohydrodynamics) model presented in Nabert et al. (2013).This model uses series expansions of density, velocity, gas pressure, and magnetic field up to the first order away from the stagnation streamline to simplify the stationary MHD equations.Consequently, the model is only valid close to the stagnation streamline.The x axis of the coordinate system is chosen to be along the stagnation streamline, the z axis is along the Earth's dipole axis, which is assumed to be orthogonal to the x axis, and the y axis completes a right-hand orthogonal system.Note that the y and z directions are called tangential hereafter.The series expansion procedure provides the following set of ordinary differential equations (Nabert et al., 2013): Here, ρ 0 represents the density, u x0 the velocity component in x direction, B z0 the magnetic field component along the Earth dipole moment, and p 0 the gas pressure.The index 0 indicates the zeroth order approximation used.The derivatives marked by the prime are with respect to the stagnation streamline direction x.The vacuum permeability is µ 0 = 4π × 10 −7 N/A 2 .The tangential velocity derivatives u y10 and u z01 , as well as the magnetic field derivative B x01 , are given by u y10 = 0.8 u/x BS , u z01 = u/x BS , and B x01 ≈ 0 (Nabert et al., 2013).Here, u denotes the x component of the velocity difference across the bow shock, which is calculated using Rankine-Hugoniot relations.The bowshock distance to the Earth's center is denoted by x BS .Figure 1 sketches the workflow of this model.
The solar wind parameters are used as input parameters to obtain the magnetosheath solution.Note that this evaluation procedure is called forward model.In a later step, yet to be discussed, these boundary values are modified to match any observed magnetosheath situation to our calculation results using a so-called reverse model approach.The solar wind parameters are the density ρ SW , velocity's x component u SW , magnetic field's z component B SW , and pressure p SW .Note that an x component of the magnetic field does not affect the solution of the zeroth order model used.This limitation of the magnetosheath model restricts our solar wind reconstruction to the y and z components of the solar wind magnetic field.The x component cannot be reconstructed using this magnetosheath model.A y component of the magnetic field can be taken into account in this model by replacing the z component B z0 by the magnitude of the tangential magnetic field, i.e., B z0 ← sgn(B z0 ) • (B 2 y0 + B 2 z0 ) 0.5 .Here, sgn() denotes the signum function.Then, the solar wind magnetic field B SW is replaced by B SW ← sgn(B z,SW ) • (B 2 y,SW + B 2 z,SW ) 0.5 .At the bow shock (x = x BS ), the solar wind values are transformed into post-shock magnetosheath conditions using Rankine-Hugoniot relations (e.g., Petrinec and Russell, 1997).The post-shock values ρ 0 (x = x BS ), u x0 (x = x BS ), B z0 (x = x BS ), and p 0 (x = x BS ) are the boundary conditions for the ordinary differential Eqs.(1)-(4) presented above.Furthermore, the bow-shock distance needs to be known in order to calculate u y10 and u z01 .The bow-shock distance x BS is the sum of the magnetopause distance to the Earth's center x MP and the magnetosheath thickness x MS .An initial estimator of the magnetopause distance denoted by x 0 MP is given  Nabert et al. (2013).The solar wind conditions are the input to the model.First, the bow-shock distance is initially estimated.Then, the Rankine-Hugoniot relations are solved and determine the postshock values, which are used as boundary conditions for the system of differential equations (DEQ).Solving this system leads to a new estimate of bow-shock distance.After some iterations, the bow-shock distance converges and the magnetosheath solution is obtained.
by Pudovkin and Semenov (1977): where the parameter of magnetopause geometry f = 2.44, the flow deflection parameter K = 0.89, and a magnetic dipole moment M = 8 × 10 15 Tm 3 are valid for a terrestrial situation.The initial magnetosheath thickness x 0 MS is estimated by an analytically derived formula (Nabert et al., 2013): where the deceleration at the shock is g u := u SW /u x0 (x = x BS ) and m BS is a measure for the solar wind magnetization defined by Here, p mag (x = x BS ) := B 2 z0 (x = x BS )/2µ 0 is the postshock magnetic pressure and γ = 5/3 is the ratio of specific heats.The sum of Eqs. ( 5) and ( 6) estimates the initial bowshock distance: Next, the system of differential Eqs. ( 1)-( 4) is solved numerically.Therefore, the differential equations are written as a matrix-vector equation: This equation can be represented by with q := (ρ 0 , u x0 , B z0 , p 0 ) T , the coefficient matrix A, and the inhomogeneity vector b.Note that b depends not only on q but also on u y10 and u z01 , initially determined by x 0 BS .For numerical calculation, the magnetosheath solution q is discretized along the x direction: where the position index n is space discretization.Consequently, Eq. (10) transforms into where x is the step size along the x direction.Equation ( 12) is solved with respect to q n+1 using a Gaussian elimination algorithm.Starting from the bow shock (x = x BS ), the solution is evaluated earthward until the magnetopause is reached, i.e., until the flow velocity vanishes (u x0 = 0).The distance in space between bow shock and magnetopause yields a new estimator of the magnetosheath thickness x 1 MS .The planetary magnetic field is approximated by a dipole with its moment along the z direction.The dipole is included via inner boundary conditions; i.e., the calculated total pressure p tot (sum of gas pressure, magnetic pressure, and dynamic pressure) at the magnetopause is related to the planetary field as follows (Nabert et al., 2013): Note that due to the iteration procedure used, the actual stagnation pressure at the magnetopause will differ from the stagnation pressure Kρ SW u 2 SW in Eq. ( 5) when initiating the iteration.Now, Eq. ( 13) is solved with respect to x MP .This gives a new estimator x 1 MP for the magnetopause distance.Together with the calculated magnetosheath thickness x 1 MS , a new estimator for the bow-shock distance is obtained: As depicted in Fig. 1, the newly estimated bow-shock distance is used to run the scheme again until the bow-shock distance converges.Then, the magnetosheath solution q n is obtained.
The magnetosheath solution is determined along the stagnation streamline.According to Nabert et al. (2013), extending the solution to locations off the stagnation streamline requires one to shift the corresponding x coordinate with respect to the paraboloid coordinates used.To transfer the spacecraft's location (x SC , y SC , z SC ) into the parabolicshifted coordinate system of the model indicated by x SC,shift , the following equation is used: where c y := c MP,y − c BS,y and c z := c MP,z − c BS,z .The curvature parameters are determined to be c BS,y = 0.4/x BS , c BS,z = 0.5/x BS , c MP,y = 0.4/x MP , and c MP,z = 0.5/x MP .However, due to the assumptions made to derive the magnetosheath model, the solution is only valid in the vicinity of the stagnation streamline.

Data assimilation
The spacecraft at a position (x SC , y SC , z SC ) in the magnetosheath measures the density ρ SC , velocity u SC , magnetic field B SC , and pressure p SC .Values for these quantities at a certain position form a single data point.Each data point is related to certain solar wind parameters individually by the model introduced above and the succession of steadystate states approximates any dynamical behavior of the magnetosheath.To find the appropriate solar wind parameters, the magnetosheath solution needs to match the data point at its position.Starting from an arbitrary first guess of the solar wind parameters, the parameters are modified until the solution matches the data, which is called data assimilation.It seems natural to start from typical solar wind conditions at Earth's orbit, such as ρ SW,t = 8.0•m P cm −3 , u SW,t = 400 km s −1 , B SW,t = 5 nT, and p SW,t = 1 × 10 −11 Pa.Note that the proton mass is m P = 1.67262178×10 −27 kg.To simplify the notation, we set m P = 1 in the following, so that the mass density is presented as a particle density.
During the process of varying the solar wind parameters, it is important to quantify the agreement between data and model.Therefore, a cost function is introduced.Consider, for example, the density ρ SC measured at (x SC , y SC , z SC ) and assume the model proposes a density ρ.The difference of these two values is a measure of the density agreement of model and data.The square of this difference is taken to obtain a function with a minimum: (ρ − ρ SC ) 2 .In order to sum up such terms for different quantities, these terms are normalized to avoid effects from different units.The measured values themselves can be taken as normalization.The cost function K sums up all contributions from the different quantities: The cost function K is minimized with respect to the solar wind parameters ρ SW , u SW , B SW , and p SW .Thereby, the step size is proportional to the value of the cost function.If the model used is able to represent the physics of the situation considered and data errors are negligible, the cost function vanishes in its (global) minimum.Taking data errors and limitations of the model into account, the global minimum of the cost function gives the best estimator for the solar wind parameters.
Starting from the typical solar wind conditions presented above, a gradient-based minimization algorithm is used to minimize the cost function.We use a steepest descent method without line search as minimization algorithm.If the gradient vanishes, a minimum of the cost function is determined.Using an adjoint approach the gradient required for this minimization scheme is determined as illustrated in the next section.

The reverse magnetosheath model -adjoint approach
An adjoint approach offers an efficient way to calculate the gradient of a cost function (e.g., Meader and Martins, 2012).
The numerical implementation of the magnetosheath model provides a source code, which calculates the cost function from the initial solar wind parameters.Here, the adjoint approach is performed with an automatic differentiation tool.Such a tool is able to transfer the source code into a new modified source code, which automatically calculates the derivative of the cost function (Wengert, 1964).Therefore, the automatic differentiation tool dissembles the code into elementary operations, i.e., simple arithmetic operations.To illustrate this, consider a simple example of a cost function similar to Rall and Corliss (1996) defined by which depends on the two parameters ρ SW and u SW .The example discussed here illustrates how the tool is processing the code.An appropriate source code implementation to calculate this function with elementary operations only is Using the chain rule, the derivative of the cost function (Eq.16) can be expressed as follows: Now, the automatic differentiation tool inserts calculation of the matrices ∂(l 1 , l 2 )/∂(ρ SW , u SW ), ∂(l 3 , l 4 )/∂(l 1 , l 2 ), and ∂l 5 /∂(l 3 , l 4 ) needed to calculate the derivative.This gives Here, the arrow to the left (e.g., a ← a • b) means that the expression on the right side (a •b) is calculated first, and then, the variable on the left side (a) is replaced by the result of the right side.The entities of the matrices are elementary derivatives such as sin(l 2 ) , which are substituted from a library used by the automatic differentiation tool: sin(l 2 ) → cos(l 2 ).This procedure has the advantage whereby the calculated derivatives are not subjected to errors of finite difference approximations because analytical expressions are inserted.
The new code (Eq.19) calculates the derivative parallel to the execution of the forward model.This procedure is called forward differentiation.However, this is not the most efficient way to calculate the derivative.It is more efficient to use the automatic differentiation in reverse mode, which is nothing else but the adjoint approach.A short introduction of the mathematical motivation of adjoint methods is given in Appendix A. Executing the code above performs the calculation of the derivative as shown in Eq. ( 18) from the right to the left.During this procedure, 12 elementary multiplications are performed.However, the calculation from the left to the right in Eq. ( 18) needs only 8 elementary multiplications.Consequently, the reverse mode is more time efficient.The transpose of Eq. ( 18) reverses the order similar to the calculation described above.Note that adjoint and transpose are equivalent in real space and this is why the reverse approach is called the adjoint approach.
The adjoint approach of the more complex cost function with the magnetosheath model presented in the previous section is derived by the OpenAD/F (open-source tool for the automatic generation of adjoint code from Fortran 95 source code) tool (Utke et al., 2008).The code contains a Fortran subroutine evaluating the magnetosheath model as shown in Fig. 1 with the solar wind parameters and the cost function as input and output variables.To prepare the subroutine for the OpenAD/F tool, solar wind parameters, such as rhoSW, are declared as independent variables by a !$openadINDE-PENDENT(rhoSW) statement.The cost function (cost) is declared as dependent by a !$openadDEPENDENT(cost) statement.The tool transfers an ordinary variable y into a structure, which contains the value of the variable y%v and the value of the derivative of the cost function with respect to y denoted by y%d.In order to enable the code to operate in reverse mode, the derivative of the cost function is initialized by 1 because the derivative of the cost function with respect to the cost function gives 1.Then, the OpenAD/F tool can be applied in reverse mode to the source code.First, the tool executes a lexical, syntactic, and semantic analysis, which results in an intermediate representation, called whirl.Then, the OpenAnalysis module of the OpenAD/F tool produces call and flow graphs.This information is combined with the whirl code into a new representation of the code called xaif.This code representation is transferred into a code that calculates the derivative parallel to the execution via the xaif-Booster module.Afterwards, the tool transfers this representation back to ordinary Fortran code.A call of the new subroutine is able to calculate the gradient.
For the four independent parameters that we consider, i.e., ρ SW , u SW , B SW , and p SW , the adjoint approach can be 4 times faster.To calculate in reverse mode requires one to execute the code before the derivatives are calculated because the derivative terms require the results calculated later in the code.Therefore, all variables used during the execution of the code to calculate the cost function need to be stored.The implementation of checkpointing, i.e., splitting the adjoint calculation into separate parts to reduce memory requirements, was not necessary for our adjoint code.

A technical problem and its solution
The cost function has a lot of very small local minima on top of the monotonous cost function with its global minimum.These minima result from the iterative modification of the bow shock as described in Fig. 1.The bow-shock distance is used to calculate u y10 and u z01 .Small modifications change the value of the cost function differently compared to the global shape of the cost function.If a finite difference approach is applied, the local minima can be ignored using a step size for the finite difference, which is larger than the extent of such a local minimum.However, it is difficult to determine the correct finite difference size because the size varies depending on the solar wind conditions.
The situation is different for an automatic differentiation approach because it is based on analytical expressions.Unfortunately, without code modifications, the gradient is calculated with respect to the local minima.Consider, e.g., the one-dimensional cost function h(x) := h q (x) + h o (x) with the quadratic function h q (x) := x 2 and the oscillation function h o (x) := (cos(10x) − 1) 2 .The global minimum is determined by the quadratic function, which decreases monotonically on both sides towards the minimum.However, the oscillation function superposes local minima to the function.The derivative ∂ x h(x) = 2x − 20 sin(10x)(cos(10x) − 1) oscillates around the global minimum due to the latter contribution ∂ x h o (x).Following the gradient usually gives a local minimum.This problem can be solved by replacing h o (x) by a sum of step functions h o (x) = i a i θ (x − x i ) with the coefficients a i and x i , which are determined to fit h o (x).Here, θ (x) denotes a step function.The analytical derivative of h o vanishes if x = x i .Consequently, ∂ x h(x) = ∂ x h q (x)+∂ x h o (x) = ∂ x h q (x) = 2x and the global minimum is found by following the gradient.In our magnetosheath model, the modifications of u y10 and u z01 produce the local minima and we excluded them by introducing step functions.This was done by rounding the tangential velocities.Note that the rounding error is chosen to be very small on the order of the numerical error, essentially not effecting the result of the cost function.Consequently, the calculation of u y10 and u z01 does not contribute to the calculated gradient and a gradient with respect to the global shape is obtained.
Figure 2 shows the resulting global shape of the cost function depending on the solar wind density ρ SW and solar wind velocity u SW .
The solar wind magnetic field and the solar wind pressure are zero.The spacecraft values are set to ρ SC = 20.4cm −3 , u SC = 48.0km s −1 , B SC = 0.1 nT, and p SC = 0.59 nPa at the spacecraft's position x SC = 12.8 R E , y SC = 0.9 R E , and z SC = 3.1 R E , with the Earth's radius R E = 6371 km.The red areas on the lower left and upper right corner are cost function values where the spacecraft's position is outside the calculated magnetosheath.Consequently, the cost function is high due to the bad matching of data and model.Note that the magnetospheric values, i.e., the values outside the magnetosheath on the earthward side, are set to the magnetopause values except from the magnetic field, which increases earthward.The middle area of the plot offers a smooth global minimum as sketched in Fig. 2, which can be determined by a gradient minimization.
The rounding procedure introduced above is used to neglect the small local minima due to iterative bow-shock modifications for the calculation of the gradient.To ensure that the minimum found by the gradient-based minimization is the local minimum, which corresponds to the global minimum, we calculate the cost function value on a parameter grid a couple of steps in each direction around the minimum.We choose the grid step sizes ρ = 0.035 cm −3 , u = 2 km s −1 , B = 0.035nT, and p = 1.5 × 10 −3 nPa.Note that the step sizes determine the accuracy of the reconstructed parameters.

Application
In the Nabert et al. (2013) magnetosheath model, the time resolution is limited because of the steady-state approximation in the model.A disturbance due to a change of solar wind conditions needs to pass the magnetosheath to obtain a stationary state.A typical subsolar magnetosheath thickness is about 2.5 R E and the average flow velocity in the magnetosheath is about 50 km s −1 .Then, the transit time for the disturbance is approximately 5 min.Consequently, we use spacecraft data with a time resolution of 5 min for our reconstruction.
The temperature is often very difficult to measure precisely.Therefore, we assume a cold plasma approximation for the solar wind.A variation of the solar wind temperature does not affect the reconstruction results much for the results presented here.Note that the temperature T 0 is related to the pressure used in the magnetosheath model via ρ 0 k B T 0 = p 0 , with the Boltzmann constant k B = 1.3806488 × 10 −23 J/K.
We have chosen THEMIS magnetosheath transitions where a second THEMIS spacecraft observes the solar wind conditions to validate our method.

Reconstruction close to the stagnation streamline
The reconstruction method is applied to magnetosheath transitions along the stagnation streamline.First, the magne- tosheath transition of the THEMIS spacecraft THC on 1 August 2009 is discussed.Starting at the magnetopause close to the stagnation point, the spacecraft traverses the magnetosheath and crosses the subsolar point at the bow shock about 2.5 h later.
The THC magnetosheath data and the best model prediction, which could be found in the magnetosheath by varying the solar wind boundary parameters, are shown in Fig. 3.Note that the data are presented in the same coordinates as the model with the z direction along the dipole axis (see Sect. 2.1).
It is seen that the observation and model prediction agree very well.Only at the magnetopause and at the bow shock do the model and the observations significantly differ.Our magnetosheath model represents the bow shock by a step function profile using Rankine-Hugoniot relations to calculate the jump conditions.This approximation is not correct close to the shock where instead a smooth transition is observed.Our algorithm fits the values within this transition region with values of either pre-or post-shock values, and thus the reconstruction fails there.The magnetopause is approximated as a rigid boundary by a vanishing flow velocity in the model.This neglects, e.g., single particle processes where high energy particles can penetrate the magnetopause boundary and diffusion processes.Our magnetosheath model is extended into the magnetosphere using the magnetopause values of velocity, density, and magnetic field.Therefore, our approach is able to reproduce the magnetospheric values observed approximately.
The results of the solar wind reconstruction are displayed together with OMNI solar wind data and the solar wind observations by the THEMIS spacecraft THB in Fig. 4.
The calculated solar wind conditions show overall a good agreement with the solar wind data.The first 10 min interval corresponds to a reconstruction from magnetospheric observations.Due to the magnetospheric extension of our magnetosheath model, the magnetosheath model is able to estimate the solar wind parameters roughly during this time interval.For example, a different choice of the solar wind magnetic field parameter would lead to a different magnetopause magnetic field, and consequently the magnetospheric field of the model does not fit the magnetospheric observations.A higher solar wind dynamic pressure is related to a different magnetopause location, and again, the model predictions would differ from the observations.However, this reconstruction method is valid only in the vicinity of the magnetopause.Our approach cannot be used to reconstruct accurate solar wind parameters far in the magnetosphere.
From 03:20 to 04:00 UT the density of the solar wind data from OMNI, THB, and the reconstructed solar wind parameters of THC differ slightly.This might be due to spacecraft potentials (e.g., McFadden et al., 2008) or due to a varying spatial density distribution of the solar wind.Note that the THB spacecraft is located about 12 R E from THC in y direction.The density variation between 05:15 and 5:40 UT ob- served by the solar wind monitor data is well reconstructed.The average value of the solar wind velocity is well reconstructed; however, the velocity reconstruction shows more variability than the solar wind data from THB and OMNI.A slight increase of the solar wind magnetic field towards the bow shock is seen in all solar wind data in the lower panel.
We ensure that our algorithm produces the same results for different initial solar wind parameters.As an example, the minimization of the cost function of the data point at 05:30 UT in Fig. 3 for three different initial values of the density is presented in Fig. 5.
A different initial choice leads to the same reconstructed solar wind parameters, only the number of iteration steps differ.Note that we have chosen a cold plasma approximation to fix the solar wind pressure.
To examine the reconstruction method in more detail, we present a second magnetosheath transition close to the stagnation streamline, which contains more solar wind variations.
On 28 September 2008, the THC spacecraft crosses the bow shock and traverses the magnetosheath.The corresponding THC observations and the model fit are shown in Fig. 6 with the bow shock on the left-hand side and the magnetopause on the right-hand side.Both, model and observations, fit very well, except for the bow-shock region similar to the previous example.The reconstruction of the solar wind parameters are compared to THB and OMNI solar wind data in Fig. 7.
One expects that the variations of the THB solar wind data are in general more precise than OMNI solar wind data because THB is located close to the subsolar bow shock and observes the local solar wind conditions.Note that THB enters the magnetosheath at 20:15 UT; therefore, afterwards no THB solar wind observations are available.The solar wind density varies between 4 and 7 cm −3 and the recon-  structed solar wind data reproduce the solar wind variations seen on timescales of about 15 min.For example, the den-sity dip around 16:25 UT, the increase and decrease around 17:55 UT, and the increase at 18:15 UT agree very well.The reconstructed density shows a small offset of about 0.2 cm −3 compared to the solar wind monitor data.Such an offset can be produced by a small spacecraft potential (McFadden et al., 2008).Comparing THC and THB solar wind observations, typical offsets of about 10 % can be observed.Thus, we assume that 10 % is a typical data error for an offset.The solar wind reconstruction fails at the bow shock and at the magnetopause as discussed before.The velocity reconstruction reproduces the solar wind variations, e.g., observed around 16:50 UT.However, similar to the previous example, there are additional small variations on timescales less than 10 min in the reconstruction, which are less pronounced or not present at all in the solar wind data, e.g., at 18:10 UT.Both solar wind data, as well as the reconstruction, provide a slightly different offset.The magnetic field data are directed approximately along the z direction during the transition and varies between 3 and 6 nT.Similar to the density variations, the magnetic field variations on scales larger than 15 min are well reconstructed from the THC magnetosheath data, e.g., the magnetic field decreases at 15:40, around 18:00, and at 19:00 UT.
Finally, a third THC magnetosheath crossing is used to investigate our reconstruction method.The THC data from 16 September 2009 were used to reconstruct the solar wind data during the magnetosheath transition.The magnetosheath data and model predictions agree very well, similar to the previous examples discussed.THB and OMNI solar wind data are compared to the THC reconstruction in Fig. 8.
The solar wind densities show a small offset.It is seen that a density and magnetic field variation at 19:35 UT is observed by THB and also seen in the reconstruction, but not in the OMNI data.This indicates that this variation is a local solar wind phenomenon.The reconstruction can reproduce such local solar wind variations better than the OMNI solar wind predictions for the subsolar point.Closer to the bow shock and magnetopause, the reconstruction fails as discussed previously.

Reconstruction alongside the stagnation streamline
The magnetosheath model used is limited to the stagnation streamline region.To evaluate the applicability of the reconstruction using observations at some distance off the stagnation streamline, we consider a magnetosheath traversal on 24 August 2008 where the THC spacecraft is positioned 5-6 R E off the stagnation streamline.The corresponding density, velocity, and magnetic field observations are shown in Fig. 9 together with the best magnetosheath model predictions for any choice of the solar wind parameters.
The density differs closer to the magnetopause location, seen on the right-hand side of the plot.The reconstruction method takes not only the density, velocity, and magnetic field magnetosheath data into account, but also the space-  craft's position where the data point was measured.Thus, more measured parameters are used (x SC , y SC , z SC , ρ SC , u SC , B SC , and p SC ) than free solar wind boundary parameters (ρ SW , u SW , B SW , and p SW ), which need to be determined.The misfit of model prediction and data can be due to limitations of the magnetosheath model itself or due to data errors.
The reconstructed solar wind parameters and the OMNI, as well as the THB solar wind data, are presented in Fig. 10.
The solar wind density of THB and OMNI show an offset of about 15 %.Because the electron and ion densities are slightly different, a spacecraft potential might be present, which can explain such an offset.The reconstructed solar wind density agrees with the THB density.However, a significantly different density behavior can be observed at 00:50 UT, where the density increases in the reconstruction but not in the solar wind monitor data.The corresponding density peak can be observed in the magnetosheath data in Fig. 9.This might be either a very local solar wind density enhancement or a time-dependent process, which is not covered by our model.Apart from that, in the first 1.5 h the differences in the reconstruction of the solar wind data are similar to the previous results for the stagnation streamline events.Closer to the magnetopause, the calculated solar wind parameters predict a continuous velocity increase during the transition.However, the OMNI and THB solar wind data offer a nearly constant solar wind velocity.
More magnetosheath transitions in the neighborhood of the stagnation streamline were investigated.Nearly all solar wind reconstructions from magnetosheath data alongside the stagnation streamline offer a systematic velocity increase.The increase is more significant far away from the stagna-tion streamline.The solar wind velocity distribution is not correlated with a magnetosheath transition of a spacecraft.Consequently, one would not expect an increase of the velocity in all transitions alongside the stagnation streamline.Thus, the systematic velocity increase is due to model limitations.Note that this limitation can be examined without solar wind monitor data because of the systematic increase in the reconstructed solar wind velocity.This is not surprising because the series expansions used to derive the magnetosheath model restricts the model to the stagnation streamline (Nabert et al., 2013).The model predicts an approximately continuous linear decrease of the velocity in the magnetosheath on the stagnation streamline.However, alongside this streamline, the velocity also decreases, but it still jumps at the magnetopause from a certain value to zero, which is outside the scope of the model used here (Nabert et al., 2013).

Conclusions
We conclude that our method using the Nabert et al. (2013) magnetosheath model is suitable to reconstruct the solar wind parameters during a magnetosheath transition.The reconstruction was validated by comparing the reconstructed solar wind parameters to solar wind spacecraft observations.The limitations of the magnetosheath model are its time resolution of about 15 min and its restriction to magnetosheath transitions close to the stagnation streamline.Note that these limitations are due to the magnetosheath model used and not due to the adjoint approach or minimization procedure.The adjoint approach presented can be applied to other spacecraft data such as CLUSTER (e.g., Escoubet et al., 2001) as well as any other magnetosheath models.For example, a global magnetohydrodynamic simulation code can be used instead.The reconstruction of solar wind parameters might be used for the data from the Messenger mission (e.g., Solomon et al., 2007) or the upcoming BepiColombo mission (e.g., Benkhoff et al., 2010) at the planet Mercury.Then, a Mercury magnetosheath model can be examined, even if in situ solar wind observations are not available.

Figure 1 .
Figure1.Scheme of the forward magnetosheath model presented inNabert et al. (2013).The solar wind conditions are the input to the model.First, the bow-shock distance is initially estimated.Then, the Rankine-Hugoniot relations are solved and determine the postshock values, which are used as boundary conditions for the system of differential equations (DEQ).Solving this system leads to a new estimate of bow-shock distance.After some iterations, the bow-shock distance converges and the magnetosheath solution is obtained.

Figure 2 .
Figure 2. Color-coded value of the cost function with respect to solar wind velocity and density.The solar wind magnetic field and gas pressure is set to zero.The global minimum can be obtained using a gradient-based minimization.

Figure 3 .
Figure 3.The magnetosheath density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel), which are observed by THEMIS on 1 August 2009, are presented in red.The best match for the data to the magnetosheath model for any initial choice of the solar wind boundary parameters was found for the magnetosheath model predictions, which are shown in blue.

Figure 4 .
Figure 4.The reconstructed solar wind boundary conditions of the model density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel) for the magnetosheath data presented in Fig. 3 are shown in blue.Additionally, the THB (red) and OMNI (green) solar wind data during this time interval is shown.

Figure 5 .
Figure 5.The value of the cost function at each iteration step is shown for three different initial densities.Independent of the starting value, the cost function converges to the identical value.

Figure 6 .
Figure 6.The magnetosheath density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel), which are observed by THEMIS on 28 September 2008, are presented in red.The best match for the data to the magnetosheath model for any initial choice of the solar wind boundary parameters was found for the magnetosheath model predictions, which are shown in blue.

Figure 7 .
Figure 7.The reconstructed solar wind boundary conditions of the model density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel) for the magnetosheath data presented in Fig. 6 are shown in blue.Additionally, the THB (red) and OMNI (green) solar wind data during this time interval are shown.

Figure 8 .
Figure 8.The reconstructed solar wind boundary conditions of the model density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel) for the magnetosheath transition on 16 September 2008 are shown in blue.Additionally, the THB (red) and OMNI (green) solar wind data during this time interval are shown.

Figure 9 .
Figure 9.The magnetosheath density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel), which are observed by THEMIS on 24 August 2008, are presented in red.The best match for the data to the magnetosheath model for any choice of the solar wind boundary parameters was found for the magnetosheath model predictions, which are shown in blue.

Figure 10 .
Figure 10.The reconstructed solar wind boundary conditions of the model density (top panel), the ion velocity (middle panel), and the magnetic field (bottom panel) for the magnetosheath data presented in Fig. 9 are shown in blue.Additionally, the THB (red) and OMNI (green) solar wind data during this time interval are shown.